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NOMENCLATURE 


a 

SS 

Sound  speed 

A, 

A 

A 

= 

Fractional  area  of  cell  side 

B, 

B 

= 

Representation  of  general  flow  quantity 

E 

= 

Specific  total  energy 
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= 

Fraction  of  cell  volume 

I, 
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Specific  internal  energy 
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r  UN, 

JB 

= 

Grid  locations  of  muzzle  exit 

IMAX, 

JMAX 

= 

Maximum  radial,  axial  grid  lines 

i  , 

j. 

k 

= 

Finite  difference  cell  indices 

P 

= 

Pressure 

M/ 

hp 

= 

Increment  in  mass  across  cell  sides, 
projectile  mass 

A 

n 

Unit  normal  vector 

N 

= 

Normal  coordinate 

r. 

Z 

= 

Radial  and  axial  coordinates 

S 

= 

Surface  area 

T 

- 

Tangential  coordinate 

Ur 

Uf 

u. 

U,  U 

= 

Radial  velocities 

V, 

•*w 

V, 

Vr 

V 

= 

Axial  velocities 

V 

= 

Volume 

VP 

= 

Projectile  axial  velocity 

e 

s= 

Covolume 

a 

= 

Increment  in  length,  time,  etc. 

xi 

Truncation  error  or  unit  vector 
Ratio  of  specific  heats 
Density 

Angle  of  normal  to  solid  boundary 


Refer  to  finite  difference  cell 
Refers  to  muzzle  exit  flow  properties 
Refers  to  projectile 

Refers  to  beads  defining  solid  boundary 

Refers  to  time  level  tn 


1.0  INTRODUCTION 


In  this  report,  a  numerical  technique  is  described 
for  simulating  the  time  dependent,  axisymmetric  flow  of 
muzzle  gases  exiting  a  gun  barrel.  The  transient  projec¬ 
tile  boundary  is  included  in  the  flowfield  in  a  dynami¬ 
cally  coupled  manner  and  the  attendant  motion  of  ambient 
air,  due  to  projectile  motion  and  muzzle  jet  expansion, 
is  described.  The  present  technique  is  based  on  a  gener¬ 
alization  of  the  fluid-in-cell  differencing  technique  of 
Rich1  as  incorporated  in  the  OIL  code  developed  by  Johnson  . 
The  modifications  include  a  method  for  treating  arbitrary 
fixed  and  moving  solid  boundaries  within  the  Eulerian  fi¬ 
nite  difference  grid.  The  method  is  developed  in  general 
for  application  to  solid  boundaries  moving  or  deforming  in 
an  arbitrary  manner  but  is  specifically  applied  in  this 
study  to  treat  the  rigid  body  motion  of  a  projectile  exit¬ 
ing  a  gun  barrel.  This  is  consistent  with  the  goal  of  the 
study,  which  was  to  develop  a  computational  tool  for  those 
interested  in  studying  the  muzzle  gas  flow  through  complex 
muzzle  devices.  The  computer  program  generated  during 
the  study  is  called  SAMS  (Small  Arms  Muzzle  Simulator)  and 
it  is  described  in  detail  in  a  separate  users  manual3. 

Transitional  ballistics  is  concerned  with  the  flow 
of  muzzle  gases  through  a  muzzle  device  and  the  interac¬ 
tion  of  that  flow  with  the  projectile.  The  description 
of  this  flow  regime  has  become  more  important  with  the 
trend  in  modern  gun  design  to  higher  muzzle  velocities, 
more  powerful  shell  loadings  and  shorter  barrels.  The 
muzzle  device,  needed  to  reduce  recoil,  flash,  noise  or 
projectile  dispersion,  has  become  an  integral  part  of  wea¬ 
pon  design.  The  complete  description  of  muzzle  gas  flow 
through  a  muzzle  device  of  complex  geometry  requires  the 
solution  of  the  partial  differential  equations  of  unsteady 
gas  dynamics  with  transient  boundaries.  The  flowfield  is 
complicated  by  a  wealth  of  fluid  mechanical  phenomena  in¬ 
cluding  shocks,  contact  surfaces,  jets  and  nonequilibrium 
flow.  Because  of  this  inherent  complexity,  numerical  solu¬ 
tion  techniques  are  clearly  indicated. 

Many  finite  difference  techniques  have  been  developed 
for  the  solution  of  fluid  dynamics  problems4'5,  which  use 
different  approaches  to  the  finite  difference  approxima¬ 
tion  of  the  fluid  conservation  equations.  The  techniques 
are  all  characterized  by  forward  differencing  in  time  and 
are  loosely  termed  time  dependent  methods  and  nearly  all 


fall  into  one  of  two  basic  categories,  depending  upon 
whether  they  are  written  in  terms  of  Lagrangian  or 
Eulerian  coordinates.  In  Lagrangian  methods,  the  finite 
difference  grid  moves  with  the  fluid,  becoming  distorted 
as  the  unsteady  flowfield  develops.  These  methods  have 
been  successfully  applied  to  fluid  mechanics  problems  in 
which  the  fluid  does  not  become  excessively  distorted 
such  as  the  interior  ballistics  problem  and  shock  tube 
problems.  The  advantage  of  the  approach  is  that  fluid 
interfaces  and  contact  surfaces  can  be  accurately  defined. 
The  great  disadvantage  is  that  if  the  fluid  is  distorted 
such  as  in  jets,  shear  layers  or  vortices,  the  finite 
difference  grid  requires  almost  continuous  rezoning  to 
keep  it  from  literally  tying  itself  into  knots.  In 
Eulerian  methods,  the  finite  difference  grid  is  fixed  in 
the  solution  space,  and  the  fluid  essentially  flows  through 
the  grid.  The  advantage  of  this  approach  is  that  all  man¬ 
ner  of  fluid  mechanical  phenomena  can  be  treated  without 
difficulty.  The  calculations  are  stable  through  large 
deformations  and  are  typically  quite  efficient  in  terms  of 
computer  time  requirements.  The  basic  disadvantage  of  this 
method  is  that  the  regions  of  the  solution  which  require 
high  density  zoning  for  accuracy  are  not  always  known 
a  priori  and  in  fact  may  change  with  time. 

The  muzzle  blast  problem,  as  formulated  in  this 
study,  has  elements  of  both  Lagrangian  and  Eulerian  char¬ 
acteristics.  For  example,  the  flow  out  of  the  muzzle 
device  is  best  treated  in  a  coordinate  system  fixed  to 
the  muzzle  whereas  a  detailed  description  of  the  flow 
around  the  projectile  requires  coordinates  moving  with 
the  projectile.  Thus  both  Lagrangian  and  Eulerian  ap¬ 
proaches  would  require  some  degree  of  modification  for 
application  to  the  muzzle  blast  problem.  A  recent  de¬ 
velopment  in  numerical  techniques,  which  is  worth  men¬ 
tioning  in  this  regard,  combines  the  features  of  both 
techniques  and  is  known  as  the  Arbitrary-Lagrangian- 
Eulerian  (ALE)  technique6 •  7 .  In  this  method,  the  finite 
difference  mesh  may  move  with  the  fluid,  be  held  fixed, 
or  move  in  any  other  prescribed  way.  The  value  of  such 
a  capability  for  application  to  the  muzzle  blast  problem 
with  both  fixed  and  moving  boundaries  is  clear. 

Although  each  of  the  above  techniques  could  in 
principle  be  applied  to  the  muzzle  blast  problem,  the 
approach  used  in  this  study  is  to  modify  an  Eulerian  tech¬ 
nique  to  account  for  the  motion  of  solid  boundaries  within 
the  finite  difference  grid.  The  basic  technique  is  em¬ 
bodied  in  OIL  code  which  can  be  described  as  a  continuous 


fluid  Eulerian  code.  It  is  written  in  fixed  Eulerian 
(rectangular  or  axisymmetric)  coordinates  and  uses  a 
fluid-in-cell  differencing  technique.  This  technique 
is  a  generalization  of  the  successful  particle-in-cell 
technique/  which  provides  for  continuous  fluid  transport 
across  cell  boundaries.  The  differencing  method  .is  clas¬ 
sified  as  a  control  volume  approach4  and  conserves  mass, 
momentum  and  energy  within  the  solution  space. 

In  the  present  study,  the  OIL  code  is  modified  to 
permit  the  calculation  of  flows  about  arbitrary  fixed  or 
moving  solid  boundaries.  Such  calculations  require  the 
treatment  of  ’’partial  cells"  which  result  when  general 
solid  boundaries  are  introduced  in  the  rectangular  finite 
difference  grid.  The  treatment  of  partial  cells  is  based 
on  a  method  proposed  by  Rich1  and  implemented  for  fixed 
boundaries  by  Gentry,  Martin  and  Daly8.  The  method  pre¬ 
sented  herein,  is  a  further  extension  of  the  formulation 
to  account  for  arbitrary  motion  of  the  boundary.  As  ap¬ 
plied  to  the  muzzle  blast  problem  the  method  thus  pro¬ 
vides  for  the  prescription  of  a  muzzle  boundary  of  rather 
general  shape  as  well  as  moving  projectile  boundary. 

Although  the  numerical  method  has  general  applica¬ 
tion  to  flows  about  fixed  or  moving  boundaries,  its  appli¬ 
cation  in  this  study  is  directed  at  the  muzzle  blast  prob¬ 
lem.  Thus  the  computer  program  developed  during  the  study 
is  tailored  specifically  for  this  problem.  This  is  consis¬ 
tent  with  the  goal  of  the  study  which  was  to  develop  a  com¬ 
putational  tool  that  is  relatively  easy  to  utilize  in  the 
muzzle  design  process.  Thus  particuxar  emphasis  was  placed 
on  the  ease  with  which  a  problem  could  be  set  up  and  on 
stable  and  efficient  operation. 

Details  of  the  numerical  method  and  its  application 
to  the  muzzle  blast  problem  are  given  in  the  following  sec¬ 
tions.  A  brief  description  of  the  muzzle  blast  flowfieid 
and  the  resulting  requirements  placed  on  the  numerical 
technique  is  given  in  Section  2.0.  In  Section  3.0,  a  de¬ 
tailed  description  of  the  numerical  method,  including  com¬ 
ment  on  stability  and  accuracy  considerations,  is  presented. 
In  Section  4.0  the  specifics  of  the  method  for  the  muzzle 
blast  problem  are  described  and  Section  5.0  presents  cal¬ 
culated  results.  Finally,  conclusions  and  recommendations 
are  presented  in  Section  6.0. 
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DESCRIPTION  OF  MUZZLE  BLAST  FLOWFIELD 


As  noted  in  the  introduction,  the  muzzle  blast 
flowfield  is  a  complex,  inherently  unsteady  and  multi¬ 
dimensional  flowfield  which  includes  numerous  fluid 
mechanical  phenomena.  The  problem  provides  a  severe 
test  of  any  numerical  simulation  technique,  due  to  the 
wide  range  of  flow  conditions  which  exist  as  the  tran¬ 
sient  flowfield  develops.  The  detailed  simulation  of  the 
problem  would,  however,  considerably  aid  attempts  to  alle¬ 
viate  undersirable  weapon  characteristics  such  as  recoil, 
noise,  flash  and  projectile  disperison.  Devices  to 
suppress  such  characteristics  are  termed  muzzle  devices 
and  have  been  designed  almost  exc]usively  by  trial  and 
error  experimental  programs.  By  way  of  introduction 
to  the  numerical  technique  and  test  cases  presented 
in  subsequent  sections,  some  background  on  muzzle  gas 
flow  problems  is  presented  in  this  section.  This  sec¬ 
tion  includes  a  qualitative  description  of  the  transient 
flowfield  and  some  insight  into  the  phenomena  of  interest. 


A  schematic  description  of  the  muzzle  blast  flow- 
field  is  given  in  Figure  1  above. 

The  muzzle  gases,  which  accelerate  the  projectileQ 
down  the  barrel,  consist  of  the  high  temperature  {^2500  K) 
and  pressure  (^500  atm)  products  of  propellant  combustion. 
As  the  projectile  leaves  the  muzzle  these  gases  are  per¬ 
mitted  to  vent  to  the  surrounding  air,  thas  initiating  a 
complex  sequence  of  events  which  equilibrate  the  muzzle 
gas  properties  with  the  ambient.  The  presence  of  the  high 
pressure  jet  of  muzzle  gas  is  communicated  to  the  ambient 
by  an  outward  rushing  shock  wave  which  outruns  the  projec¬ 
tile  but  decays  rapidly  in  intensity  and  velocity.  This 
weakening  "blast  wave"  is  the  source  of  much  of  the  far- 
field  sound  resulting  from  the  firing  sequence.  In  the 
meantime,  the  muzzle  gases  continue  to  expand  and  accel¬ 
erate  into  the  shock  processed  air  between  the  projectile 
and  barrel,  and  around  the  projectile.  It  is  clear  that 
while  the  projectile  is  in,  or  in  close  proximity  to, 
the  muzzle  device,  the  moving  projectile  boundary  will 
have  an  important  effect  on  the  flow  in  the  device.  In 
fact,  the  projectile  is  the  obstruction  which  causes  the 
muzzle  gas  to  expand  radially  if  permitted  by  muzzle  geom¬ 
etry.  The  jet  of  muzzle  gases  expands  around  the  projec¬ 
tile,  so  that  the  projectile  is  initially  flying  backwards 
relative  to  the  muzzle  gases.  Thus  the  projectile  is  still 
accelerating  relative  to  the  barrel,  but  the  muzzle  gases 
begin  to  decelerate  so  that  the  projectile  then  penetrates, 
in  turn,  the  gas  cloud  and  the  shock  wave  as  it  proceeds 
on  to  the  target.  The  scenario  described  above  all  hap¬ 
pens  on  a  time  scale  of  the  order  of  100  microseconds. 

Various  muzzle  devices  are  used  to  tailor  the 
muzzle  gas  flow.  Recoil  forces  can  be  significantly 
decreased  by  venting  the  muzzle  gases  radially  away  from 
the  barrel  before  the  projectile  exists.  The  tendency 
of  a  weapon  to  climb  can  be  reduced  at  the  same  time  by 
directing  the  muzzle  gases  upward.  It  is  also  believed 
that  projectile  dispersion  can  be  reduced  somewhat  by 
tailoring  the  geometry  of  the  flow  deflector.  As  noted 
above,  the  projectile  is  initially  flying  in  an  unstable 
attitude  (backwards)  through  the  muzzle  gas  cloud.  De¬ 
flecting  the  flow  away  from  the  axis  would  decrease  the 
flow  density  and  velocity  through  which  the  projectile 
is  passing,  thereby  decreasing  the  intensity  and  duration 
of  the  unstable  forces  on  the  projectile.  This  could  be 
especially  important  for  multiple  fleschette  projectiles 
which  are  considerably  less  stable  in  a  reverse  flow  than 
normal  projectiles. 


Other  undesirable  characteristics  of  the  muzzle 
flow  are  ioise  and  muzzle  flash.  A  source  of  noise  other 
than  the  initial  air  shock  wave  and  the  supersonic  pro¬ 
jectile  shock  system  ("ballistic  crack")  are  the  fluctu¬ 
ations  in  the  shear  layers  of  the  muzzle  jet  as  it  exits 
the  muzzle.  Muzzle  geometry  can  significantly  affect  the 
nature  (Intensity  and  directionality)  of  such  sources  of 
farfield  noise  by  decreasing  the  magnitude  of  flew  pres¬ 
sure  and  velocity  before  the  gas  is  vented  from  the  de¬ 
vice.  Muzzle  flash  is  caused  by  the  high  temperature, 
nonequilibrium  nature  of  the  gas  as  it  exits  the  muzzle. 
Again  muzzle  geometry  can  be  chosen  so  as  to  expand  the 
gas  to  a  temperature  below  propellant  ignition  tempera¬ 
ture  before  permitting  it  to  exit  from  the  device.  One 
method  to  accomplish  this  would  be  to  define  the  device 
geometry  in  such  a  way  that  cool  ambient  air  is  entrained 
by  the  jet  of  muzzle  gases. 

All  of  the  undesirable  weapon  characteristics  men¬ 
tioned  above  are  strongly  dependent  upon  muzzle  device 
geometry.  To  be  useful  in  any  systematic  study  of  de¬ 
vices  to  reduce  or  eliminate  such  characteristics,  a  nu¬ 
merical  simulation  technique  must  necessarily  account,  for 
complex  muzzle  geometry  and  a  moving  projectile  boundary. 
Projectile  geometry  could  also  be  important  to  the  disper¬ 
sion  caused  by  the  flow  of  muzzle  gases  over  the  projec¬ 
tile  and  could  be  studied  in  the  same  manner.  The  numer¬ 
ical  technique  and  computer  program  developed  in  the  pre¬ 
sent  study  thus  emphasize  geometry  considerations  to  fa¬ 
cilitate  such  studies. 

The  numerical  technique  for  simulating  the  complex 
flowfield  described  above  must  be  stable  through  wide  var¬ 
iations  in  flow  velocity  from  supersonic  to  stagnation 
and  in  flow  densities  from  many  times  ambient  in  the  muz¬ 
zle  jet  to  considerably  below  ambient  in  the  expansion 
following  the  air  shock.  This  coupled  with  the  need  for 
defining  complex  muzzle  geometry  indicated  that  an  Euler- 
ian  differencing  technique  would  be  most  appropriate.  The 
development  of  the  technique  used  in  this  study  is  given 
in  the  following  sections. 
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3.0 


NUMERICAL  SOLUTION  METHOD 


The  main  complexity  in  the  numerical  calculation 
of  muzzle  blast  flowfields,  described  above,  involves  the 
treatment  of  the  irregular  solid  boundaries  defining  the 
muzzle  device  and  moving  projectile.  The  importance  of 
the  computational  treatment  of  boundary  conditions  in 
any  finite  difference  technique  cannot  be  overstated,  as 
the  treatment  can  affect  not  only  the  accuracy  but  the 
stability  of  the  numerical  technique4.  Thus  the  main 
thrust  of  the  present  effort,  and  its  primary  accomplish¬ 
ment,  is  in  the  area  of  a  finite  difference  approximation 
to  irregular  fixed  and  moving  solid  boundaries.  The  nu¬ 
merical  technique  is  described  in  detail  in  this  section. 

Sections  3.1  and  3.2  below  summarize  the  relevant 
conservation  equations  and  the  fluid-in-cell  differencing 
technique  used  in  the  OIL  program.  Modifications  for 
fixed  and  moving  solid  boundaries  are  given  in  Section 
3.3  and  a  discussion  of  the  stability  and  accuracy  of  the 
overall  numerical  technique  is  presented  in  Section  3.4. 


3.1  Conservation  Equations  for  Inviscid,  Compressible 
Fluid  Mechanics  and  Solution  Methodology 


In  the  treatment  of  the  muzzle  blast  problem  in 
this  report,  the  flowfield  generated  by  the  muzzle  gases 
is  considered  within  the  context  of  inviscid,  unsteady 
gas  dynamics  of  a  single  fluid  with  real  gas  equation  of 
state.  Assuming  axisymmet* ic  flow,  the  governing  differ¬ 
ential  equations  are  written  in  cylindrical  (r,  z)  coordi¬ 
nates  as: 
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where  p,  P,  E  are  the  flxiid  density,  pressure  and  total 
energy  respectively  and  u,  v  are  the  velocity  components 
in  the  radial,  r,  and  axial,  z,  coordinate  directions  re¬ 
spectively,  and  t  is  time.  The  total  energy  of  the  fluid 
is  defined  in  the  usual  manner: 


E  =  i(u2+v2)  +  I 


(4) 


where  I  is  the  specific  internal  energy.  To  effectively 
close  the  above  system  of  equations,  a  fluid  equation  of 
state  is  appended  to  define  the  pressure  as: 


P  =  P(P,I)  (5) 


As  noted  earlier,  the  fluid-in-cell  differencing 
technique  for  approximating  the  above  equations  uses  a 
control  volume  approach,  so  that  it  is  useful  in  describ¬ 
ing  the  solution  method  to  refer  to  the  control  volume 
form  of  the  conservation  equations  which  are  written  as 
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where  V  is  an  arbitrary  volume  in  the  fluid  with  surface 
area  S  and  ri  is  the  unit  outward  normal  on  the  surface. 


In  the  finite  difference  approximation  to  the  conserva¬ 
tion  equations,  presented  in  the  next  section,  the  inte¬ 
gration  volumes  are  the  cells  of  the  finite  difference 
mesh.  The  control  volume  equations  are  thus  helpful  in 
understanding  the  treatment  of  the  transport  or  convec¬ 
tive  terms  in  the  numerical  technique. 

The  method  for  solving  the  above  system  of  equa¬ 
tions  is  briefly  summarized  here  by  way  of  introduction 
to  the  finite  difference  approximations  presented  in  the 
next  section.  The  calculations  necessary  to  advance  the 
solutions  one  stew  in  time,  At,  are  separated  into  two  dis¬ 
tinct  phases.  The  first  phase  consists  of  an  explicit 
Lagrangian  calculation  for  a  fluid  particle  defined  by 
a  computational  cell  in  which  fluid  properties  (p,  u,  v,  E) 
are  assumed  uniform.  Neglecting  convection  terms,  the 
axi symmetric  form  of  the  conservation  equations  are: 
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An  alternate  form  of  the  energy  equation  which  is  used  is: 


3 (ru)  j 

3r  j  (12) 


These  equations  are  used  in  phase  one  to  update  u,  and  I 
to  an  intermediate  state.  At  this  stage  in  a  Lagrangian 
technique  the  grid  points  defining  a  finite  difference 
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cell  (or  fluid  particle)  would  be  moved  to  their  advanced 
time  positions  and  new  cell  densities  would  be  calculated. 
In  an  Eulerian  method  as  used  here,  however,  the  grid 
points  remain  fixed  in  space  so  that  the  transport  terms 
must  be  accounted  for.  Thus  in  phase  two  of  the  calcula¬ 
tion,  mass,  momentum  and  energy  are  transported  across 
the  fixed,  cell  boundaries  thereby  completing  the  process 
of  updating  flow  quantities  to  time  t  +  At,  The  control 
volume  form  of  the  conservation  equations,  neglecting  pres¬ 
sure  terms  which  were  accounted  for  in  phase  one,  is  most 
appropriate  in  this  phase.  Each  cell  is  considered  as  a 
control  volume  and  transport  to  or  from  adjacent  cells  is 
accounted  for.  In  some  sense,  this  phase  is  equivalent 
to  rezoning  the  Lagrangian  cells,  or  in  other  words  mov¬ 
ing  the  Lagrangian  cell  boundaries  back  to  their  positions 
at  the  beginning  of  the  cycle.  This  is  used  to  advantage 
in  the  arbitrary  Lagrangian-Eulerian  technique,  in  which 
cell  boundaries  are  moved  in  a  general  manner. 

The  separation  of  the  computational  cycle  into  a 
Lagrangian  phase  and  a  convective  or  rezone  phase  origi¬ 
nated  in  the  particle-in-cell  numerical  method9,  and  has 
been  used  in  many  hydrodynamic  computer  codes,  besides 
the  OIL  code.  The  rationale  behind  the  separation  into 
phases  is  basically  that  the  equations  of  motion  are  in¬ 
herently  Lagrangian  and  are  quite  straightforward  exten¬ 
sion  of  Newtonian  dynamics  to  a  continuum.  The  technique 
thus  separates  the  physics  of  the  problem  from  the  mathe¬ 
matical  complexity  of  the  transformation  to  Eulerian  coord¬ 
inates.  In  addition  to  this  conceptual  simplification,  the 
separation  of  Lagrangian  effects  from  transport  effects 
has  the  important  benefit  of  permitting  a  finite  differ¬ 
ence  scheme  which  "globally"  conserves  mass,  momentum, 
and  energy.  For  example,  energy  is  conserved  throughout 
the  finite  difference  grid  as  a  result  of  this  separation 
into  phases  plus  the  proper  choice  of  time  centering  in 
the  Lagrangian  phase.  The  finite  difference  approximation 
to  the  conservation  equations  used  in  each  cycle  are  now 
summarized  in  the  next  section. 


3 . 2  Fluid-in-Cell  Finite  Difference  Approximation 


The  fluid-in-cell  differencing  technique,  used  in 
the  OIL  program,  Is  described  in  detail  by  Rich1,  Gentry, 
et  al.8,  and  Johnson2.  The  finite  difference  forms  used 
in  the  method  are  summarized  in  this  section  and  the 
reader  is  referred  to  the  above  references  for  a  more 
detailed  description. 
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Figure  2,  Schematic' of  Local  Cell  Structure 


In  the  finite  difference  approximation  to  the  con¬ 
tinuum  conservation  equations,  the  region  of  fluid  under 
consideration  is  covered  by  a  fixed  grid  of  rectangular 
cells  with  cell  indices  and  mesh  spacing  as  indicated  in 
Figure  2.  The  index  k  is  related  in  a  straightforward 
manner  to  i  and  j  and  denotes  the  cell  number.  Arbitrary 
cell-spacing,  in  both  r  (radius)  and  z  (axial)  directions, 
is  used  to  permit  concentration  of  cells  in  regions  of 
the  flow  in  which  a  high  degree  of  resolution  is  desired. 
For  the  purposes  of  differencing,  all  flow  quantities 
(p,  P,  I,  u,  v)  are  defined  at  cell  centers.  An  unsteady 
calculation  is  initiated  by  defining  flow  properties  in 
each  cell  at  an  initial  time  t  =  0.  The  flow  quantities 
are  then  updated  in  a  cyclic  manner  to  later  times  using 
the  finite  difference  analog  to  the  fluid  conservation 
equations.  Considering  the  fluid  state  at  some  time  t, 
the  equations  used  to  update  all  fluid  quantities  to  an 
advanced  time  t  +  At  are  now  described. 

It  is  recalled  that  the  cycle  begins  with  a  pure 
Lagrangian  calculation  (phase  one)  in  which  the  convective 
terms  in  the  Eulerian  form  of  the  equations  of  motion 
are  neglected.  The  differential  form  of  the  equations 
in  the  cylindrical  coordinates  of  interest  were  given 
in  the  previous  section  for  the  two  components  of  momen¬ 
tum  (r,  z)  (Equations  9  and  10,  respectively)  and  for 
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internal  energy  (Equation  12) .  Time  derivatives  are  ap 
proximated  in  the  usual  explicit  form,  for  example: 


(13) 


A  tilde  (M  on  all  flow  quantities  denotes  the  value  at 
the  end  of  phase  one  and  superscript  n  denotes  the  value 
at  time  t.  The  pressure  gradients  are  approximated  by  a 
leapfrog  differencing  scheme  in  which  the  pressure  on 
cell  boundaries  is  taken  as  the  simple  average  of  the 
pressure  in  the  cell  and  its  neighbor.  The  resulting  fi 
nite  difference  equations  for  updating  flow  quantities 
in  phase  one  are  written  as: 


where 
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The  use  of  the  average  velocities  u  and  v  in  the  internal 
energy  equation  is  required  for  energy  conservation  in  the 
finite  difference  analog.  Rich1  gives  a  detailed  explana¬ 
tion  for  their  use  which  will  not  be  repeated  here.  It 
is  noted  that  density  is  not  updated  in  phase  one  since 
the  mass  in  a  cell  (fluid  particle)  is  assumed  constant 
in  the  Lagrangian  sense. 
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Figure  3,  Numbering  System  for  Cell  Boundaries 
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The  density  and  other  flow  properties  are  next  up¬ 
dated  to  the  advanced  time  t  +  At  in  phase  two  which  ac¬ 
counts  for  the  transport  of  mass,  momentum  and  energy 
across  cell  boundaries,  neglected  in  phase  one.  The  num¬ 
bering  system  for  jell  boundaries  is  given  in  Figure  3  and 
used  in  the  following  discussion.  The  control  volume  form 
of  the  conservation  equations  (Equations  6,  7  and  3)  is 
applied  to  each  cell  and  flow  of  mass,  momentum  and  energy 
across  cell  boundaries  is  accounted  for.  The  finite  dif¬ 
ference  approximation  to  the  continuity  equation  is: 


amJ  +  AM^  +  AM^ 


(18) 


where  Vi,j  is  the  cell  volume  and  AM^  is  the  increment 
in  mass  which  flows  across  cell  side  i  during  the  time 
At.  The  mass  flow  across  side  2,  for  example,  is  approx¬ 
imated  using  donor  cell  differencing  by: 
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It  is  noted  that  the  density,  in  this  expression,  is  taken 
as  the  value  in  the  donor  cell  and  that  U2  is  a  weighted 
velocity  for  fluid  flowing  from  the  donor  cell.  This  ve¬ 
locity  is  determined  by  a  Taylor  series  expansion  about 
the  pell  boundary,  assuming  that  the  fluid  moves  a  distance 
A  =  U2At  in  the  time  increment  At.  The  expansion  is: 
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Solving  this  expression  for  gives : 


.'v 

(Ui,j 


'u  . 

-  u.  ,  . ) 
1-1/3 


Ari-1 


(20) 


The  use  of  this  weighted  velocity  (of  which  Equation  20 
is  just  one  of  many  approximations  tnat  could  be  used)  is 
an  essential  ingredient  to  the  donor  cell  differencing 
scheme  and  has  been  found  by  Rich1  to  have  superior  dif¬ 
fusive  properties. 

The  increments  in  mass  flow  on  the  other  sides  of 
the  cell  are  determined  in  an  equivalent  manner. 

Transport  of  momentum  and  total  energy  are  treated 
in  the  same  way  and  the  result  is: 
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where  B  denotes  u,  v  or  E.  It  is  noted  that  the  B^(i  =  1, 
2,  3,  4)  are  also  determined  by  the  respective  flow  quant¬ 
ity  in  the  donor  cell .  The  process  for  updating  all  quant¬ 
ities  to  an  advanced  time  (t  +  At)  is  completed  by  calcu¬ 
lating  the  specific  internal  energy  and  pressure  in  each 
cell  as  follows: 


The  solution  method  described  in  this  section  has 
been  determined  by  many  years  of  numerical  experimenta¬ 
tion.  Details  of  the  approximations  such  as  the  donor 
cell  (or  upstream)  differencing  and  the  velocity  weight¬ 
ing  of  mass  flux  across  cell  boundaries  are  found  to  be 
necessary  to  the  stability  and  accuracy  of  the  scheme. 

The  scheme  is  conservative  for  the  rectangular  mesh  sys-^ 
terns  used,  as  described  by  Rich1  and  Johnson2.  The  dif¬ 
ference  approximations  used  are  first  order  accurate  with 
regard  to  time  increment  and  mesh  spacing,  and  could  be 
improved  in.  a  number  of  ways.  However,  the  method  des¬ 
cribed  here  has  been  used  successfully,  and  with  good 
comparison  to  experimental  or  alternate  numerical  results, 
on  a  wide  range  of  fluid  mechanics  problems.  The  treat¬ 
ment  of  boundary  conditions,  essential  to  the  mu2zle  blast 
problem,  is  described  in  the  next  section.  This  treatment 
is  consistent  with  the  basic  scheme  discussed  here. 
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3.3 


Treatment  of  Irregular  Solid  Boundaries 


The  main  modification  to  the  OIL  code,  made  during 
this  study,  was  to  add  the  capability  of  defining  fixed 
and  moving  boundaries  of  arbitrary  shape  within  the  fi¬ 
nite  difference  solution  field.  This  modification  is  of 
prime  importance  to  the  muzzle  blast  problem  because  of  the 
professed  desire  of  predicting  the  velocity  and  pressure 
field  in  and  around  the  muzzle  device  and  projectile.  The 
treatment  of  boundary  conditions  is  generally  the  most  im¬ 
portant  and  often  the  most  difficult  aspect  of  any  numeri¬ 
cal  technique  for  fluid  mechanics  problems.  The  most  accu¬ 
rate  and  straightforward  treatment  of  solid  boundaries  is 
to  use  an  orthogonal  coordinate  system  which  matches  the 
boundary  contour4 .  This  would  be  impossible  for  the  muz¬ 
zle  blast  problem,  however,  due  to  the  existence  of  closed 
contours  (projectile  or  baffles,  etc.)  in  the  flowfield. 

The  method  described  here  is  to  overlay  the  boundary  sur¬ 
face  on  the  rectangular  grid  geometry.  A  special  account¬ 
ing  is  thus  made  of  the  "partial  cells"  on  the  boundary 
in  terms  of  both  their  geometry  and  dynamic  boundary  con¬ 
dition.  The  method  is  based  upon  the  work  of  Rich1  for 
fixed  boundaries  (Section  3.3.1)  and  is  generalized  here 
to  account  for  arbitrary  boundary  motion  in  Section  3.3.2. 


3.3.1  Fixed  Boundaries 


The  OIL  code  is  written  for  rectangular  cell  geom¬ 
etry  with  general  grid  spacing.  It  was  originally  devel¬ 
oped  for  the  calculation  of  hypervelocity  impact  and  in¬ 
tense  explosions,  so  that  no  provision  is  made  for  general 
boundaries  in  the  grid.  Later  versions  of  the  OIL  code 
account  for  fluid  interfaces  within  the  grid10'11,  per¬ 
mitting  the  calculation  of  multi-fluid  problems.  The 
treatment  used  for  such  interfaces  is  however  inappropri¬ 
ate  for  the  treatment  of  solid/fluid  boundaries  within 
the  context  of  inviscid  fluid  mechanics.  The  physical 
boundary  condition  for  inviscid  flow  (slip  boundary  condi¬ 
tion)  over  a  solid  boundary  is  that  the  normal  velocity 
of  the  fluid  at  the  boundary  equals  the  normal  velocity 
of  rhe  boundary.  For  fixed  boundaries  this  means  that 
the  flow  velocity  normal  to  the  boundary  must  equal  zero. 

The  difference  equations  given  in  the  previous  sec¬ 
tion  are  valid  only  for  interior  cells  not  on  a  boundary 
and  must  therefore  be  modified  to  account  for  the  zero 
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normal  velocity  boundary  condition.  In  the  simple  case 
in  which  the  solid  boundaries  fall  on  grid  lines  the  de¬ 
sired  boundary  condition  can  be  obtained  by  the  use  of 
fictitious  cells  inside  the  solid  boundary.  For  example, 
consider  a  cell  (i,j)  which  is  adjacent  to  a  solid  bound¬ 
ary  which  falls  on  the  i  grid  line  or  right  hand  side 
of  the  cell.  A  fictitious  cell  (i+l,j)  is  used  whose 
flow  conditions  are  completely  determined  by  the  condi¬ 
tions  in  cell  (i,j).  That  is,  in  the  normal  "reflective" 
boundary  condition  sense,  the  flov;  properties  (density, 
pressure,  internal  energy)  in  the  fictitious  cell  must  be 
equal  to  the  corresponding  flow  properties  in  cell  (i,j) 
at  all  times.  Also,  the  normal  velocity  (u  in  this  exam¬ 
ple)  in  the  cell  (i+l,j)  is  the  negative  of  the  normal 
velocity  in  cell  (i,j)  and  the  tangential  velocity  (v)  is 
at  all  times  equal  in  both  cells.  Thus  at  the  end  of 
phase  one  set: 


n 


o  =  on  P  =  p  ? 

pi+l, j  Pi, j '  Fi+l,j  Pi,j'  i+l, j  i, j 


(24) 


i+l,  j 


'Xi  ^  <v> 

.,v.t1  .=v.  . 
1,3  1+1,3  i,3 


and  at  the  end  of  phase  two  set: 


n+1  _  n+1  n+1  _  „n+l  _  ,n+l 

Pi+1, j  Pi,j'  i+l, j  Fi,j'  J'i+l,j  Ai, j 


n+1  n+1  n+1  _  „n+l 

u .  . ,  .  =  -  u .  . ,  v .  , .  .  =  v.  . 

1+1,3  1,3  i+l,3  1,3 


(25) 


Solid  boundaries  on  other  sides  of  the  grid  are  treated 
in  a  similar  manner.  If  the  solid  boundary  has  a  corner, 
as  shown  in  Figure  4  below,  the  properties  of  the  ficti¬ 
tious  cell  (i+l,j)  are  assigned  values  depending  upon 
whether  the  cell  (i , j )  or  (i+l,j+l)  is  being  considered. 
The  above  prescription  essentially  fixes  the  gradients 
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of  all  flow  quantities  equal  to  zero  at  the  surface  as 
well  as  setting  the  normal  velocity  equal  to  zero.  In 
a  mathematical  sense  this  overdetermines  the  problem.  It 
could  be  improved  by  extrapolating  flow  quantities  to  the 
boundary  based  on  more  cells  than  just  the  boundary  cell 
but  the  method  used  here  is  consistent  with  the  first  order 
differencing  scheme  and  hoped  for  increases  in  accuracy 
based  on  an  extrapolation  procedure  are  seldom  realized 
in  practice. 


Figure  4,  Schematic  of  Boundary  Cells 


The  method  described  above  for  solid  boundaries  on 
grid  lines  is  used  in  the  present  study  for  defining  the 
straight  portion  of  the  gun  barrel  "upstream"  of  the  muz¬ 
zle  device.  A  more  complicated  method  is  required  for  a 
general  body  surface  such  as  a  muzzle  device  boundary. 
Consider  a  general  curved  surface,  as  shown  in  Figure  5, 
defined  by  a  series  of  points  ("beads")  connected  by 
straight  line  segments.  The  partial  cells,  generated  by 
cuts  across  grid  lines,  are  characterized  by  five  quantities 
Fi,  j  t  (Al)  j  ,  (A2)i,j/  (A4)ifj.  Fi,i  ^e 

fractional  volume  of  the  partial  cell  and  j  (k  =  1, 

2,  3,  4)  are  the  fractional  areas  open  to  the  flow  for 
cell  sides  as  denoted  in.  Figure  5.  These  quantities 
are  defined  in  a  purely  geometrical  fashion  based  on  bead 
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location  and  the  local  grid  dimensions.  The  details  of 
the  method  for  defining  partial  cell  geometry,  for  a  gen' 
eral  fixed  or  moving  boundary,  are  given  in  Appendix  A. 

Partial 


Figure  5.  Partial  Cells  for  a  General  Curved  Surface 


The  Lagrangian  equations  used  in  phase  one  to  up¬ 
date  u,  v  and  I  are  rewritten  to  take  account  of  the  geom¬ 
etry  of  the  partial  cells  as  follows: 


Some  comment  is  required  on  the  degree  of  approxi¬ 
mation  of  the  partial  cell  geometry  inherent  in  Equations 
(26)  through  (28).  For  simplicity,  consider  just  the  axial 
momentum  equation.  Using  the  Lagrangian  form  of  the  con¬ 
trol  volume  momentum  equations  (Equation  7  with  no  trans¬ 
port  terms) ,  it  can  be  shown  that  the  axial  momentum  equa¬ 
tion  to  first  order  can  be  written  as: 
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where  Pg  is  the  fluid  pressure  along  the  solid  boundary. 
Many  approximations  to  Pg  could  be  used  which  are  accu¬ 
rate  to  the  same  order  as  the  basic  differencing  scheme. 
Equation  (27)  results  from  the  assumption  that  Ps  equals 
the  pressure  on  the  cell  side  with  minimum  opening  to  the 
flow.  This  is  the  simplest  approximation  possible  and  in 
fact  was  used  with  some  success  by  Rich1 ,  which  is  the  ba¬ 
sic  reason  for  using  it  in  this  study.  Another  possible 
approximation  would  be  to  assume  that  Ps  =  Pi,j  in  the 
normal  reflective  boundary  condition  sense.  A  potentially 
more  accurate  approximation  would  be  to  use  the  momentum 
equation  normal  to  the  solid  surface  with  a  one  sided  dif¬ 
ference  approximation  for  the  pressure  gradient  normal  to 
the  wall,  to  determine  Pg.  Unfortunately  these  alternative 
techniques  were  not  examined  in  the  present  study  but  will 
hopefully  be  the  subject  of  future  research. 

As  noted  above,  the  finite  difference  approximations 
given  here  are  just  one  of  many  alternative  forms  which  are 
consistent  with  the  first  order  accuracy  of  the  general 
difference  scheme.  The  approximation  could  be  improved  in 
a  number  of  ways,  one  of  which  would  be  to  extrapolate  flow 
quantities  into  the  fictitious  boundary  cells  as  mentioned 
earlier.  A  more  important  improvement  would  be  to  expli¬ 
citly  account  for  the  shifting  center  of  mass  of  the  par¬ 
tial  cells  in  the  Lagrangian  phase  of  the  calculation. 
Improvements  such  as  these  would  require  extensive  numer¬ 
ical  experimentation  to  verify  the  hoped  for  increase  in 
accuracy  and  were  beyond  the  scope  of  this  initial  study. 
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The  treatment  of  transport  terms  for  partial  cells 
in  phase  two  of  the  calculation  are  a  straightforward  ex¬ 
tension  of  the  general  expressions  to  account  for  the 
partial  cell  geometry.  The  partial  cell  equivalents  of 
the  mass,  momentum  and  energy  equations  (Equations  18  and 
21)  in  phase  two  are: 


3^*^  =  p?  .  +  % -  A.AM^  -i-  A0AM? 

1,3  1,3  Fi/j Vi,j  l  1  1  22 


+  a3am5  -  a4amJ 


! 


(31) 


and 


B?**  -  B?  .  +  — tt— 

1/3  1/3  0n+1F 

pi,j  i,3  1/3 


hrr  f 
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+  B„AnAM 


(32) 
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-  b4a4a<( 


where,  as  before,  B  denotes  u,  v  or  E,  V^j  is  the  rec¬ 
tangular  cell  volume  and  subscript  k  (k  =  i,  2,  3,  4)  re¬ 
fers  to  the  (i,j)  cell  sides  as  given  in  Figure  5.  The 
mass  transport,  velocity  weighting  scheme  and  donor  cell 
differencing  scheme  are  the  same  as  given  in  the  previous 
section. 

It  is  reiterated  that  the  treatment  of  fixed  bound 
aries  described  in  this  section  is  consistent  with  the 


-23 


first  order  accuracy  of  the  basic  fluid-in-cell  scheme* 
Improvements  to  the  method,  such  as  an  accurate  definition 
of  boundary  cell  center  of  mass  and  interpolation  on  the 
center  of  mass  to  determine  quantities  on  the  cell  bound¬ 
aries  will  hopefully  be  considered  in  continued  research 
on  the  method.  It  should  be  mentioned  that  the  less  exact 
method  given  here,  however,  is  at  least  as  accurate  as  the 
donor  cell  or  upwind  differencing  required  for  stability 
of  the  method.  Thus,  it  is  not  clear  that  such  improve¬ 
ments  would  lead  to  a  more  accurate  scheme. 


3.3.2  Arbitrary  Moving  Boundaries 


The  treatment  of  moving  boundaries  within  the  Eu- 
lerian  grid  is  the  essential  contribution  of  the  present 
code  development  effort.  This  permits  the  inclusion  of 
the  projectile  in  the  flowfield  which,  as  noted  earlier, 
is  important  to  the  overall  development  of  the  muzzle  gas 
flow.  The  treatment  given  here  is  a  generalization  of 
that  just  given  for  fixed  boundaries  with  provision  for 
general  time  dependent  motion  of  the  boundary.  Although 
the  technique  is  developed  for  general  motion,  it  is  spe¬ 
cialized  where  applicable  for  the  rigid  body  motion  of 
the  projectile. 

Consider  a  partial  cell,  occupied  at  a  time  t  by 
a  solid  boundary  which  is  moving  in  a  general  manner.  Tne 
boundary  motion  is  defined  by  the  motion  of  a  point  q[(Ug(t)„ 
Vq ( t ) ]  on  its  surface,  as  shown  schematically  in  Figure  6. 

As  in  the  case  of  the  fixed  boundaries,  the  moving  bound¬ 
ary  is  defined  by  straight  line  segments  joining  a  sequence 
of  beads  moving  with  its  surface.  The  position  of  the 
beads  in  the  finite  difference  grid  [rg(t),  Zq(t)]  at  any 
instant  determines  the  geometrical  factors  [Ff  j,  Ai  (i  = 

1,  2,  3,  4)]  for  the  partial  cells  all  along  the  solid 
boundary.  Again  the  details  of  the  calculation  of  these 
factors  is  given  in  Appendix  A.  The  angle  (61)  of  the 
local  normal  to  the  boundary  is  an  additional  factor,  not 
used  for  fixed  boundaries,  but  which  is  needed  for  the 
dynamics  of  moving  boundaries.  For  present  purposes,  this 
angle  is  defined  in  an  approximate  manner  as  the  angle  be¬ 
tween  the  z  axis  and  the  normal  to  the  line  segment  join¬ 
ing  the  points  at  which  the  boundary  cuts  the  cell  sides, 
as  shown  in  Figure  6. 
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where  N  and  T  are  the  normal  and  tangential  coordinates 
respectively. 

The  treatment  for  the  partial  cells  on  the  moving 
boundary  used  here  is  to  assume  that  the  normal  ve^city 
of  the  fluid  in  the  cell  is  equal  to  the  normal  velocity 
of  the  boundary  or: 


Un  sin  6!  .  +  Vn  cos  6! 
q  i,j  q  i,j 


(34) 


where  i,j  refers  to  the  partial  boundary  cell,  U„,  V„  (see 
Figure  6)  are  the  velocity  of  the  boundary  bead  at  time  t 
and  the  tilde  (^)  refers  to  values  at  the  intermediate 
state  after  phase  one.  It  is  noted  that  this  is  equiva¬ 
lent  to  prescribing  the  reflection  boundary  condition  nor¬ 
mal  to  the  surface  or: 


(35) 


Phase  one  is  completed  by  updating  the  tangential  flow 
velocity  in  the  partial  cell  based  on  a  finite  difference 
approximation  to  the  tangential  equation  of  motion.  This 
involves  approximating  the  tangential  pressure  gradient 
by  resolving  the  radial  and  axial  pressure  gradients  in¬ 
to  the  tangential  direction.  The  radial  and  axial  pres¬ 
sure  gradients  are  calculated  in  the  same  manner  as  for 
fixed  boundary  partial  cells  as  used  in  Equations  (26) 
and  (27)  above.  The  resulting  finite  difference  approxi¬ 
mation  to  the  tangential  pressure  gradient  is: 
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where  as  before 
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and  cell  sides  are  as  numbered  in  Figure  6.  Also  the 
pressures  in  fictitious  cells  completely  within  the  solid 
boundary  are  set  equal  to  the  pressure  in  the  cell  (i,j) 
as  prescribed  by  the  reflection  boundary  condition.  Based 
on  this  form  of  the  tangential  pressure  gradient,  the 
tangential  flow  velocity  is  updated  to  the  intermediate 
stage  as  follows: 


(38) 


The  axial  and  radial  velocities  i£  the  partial  cell  are 
found  simply  by  resolving  and  qT  in  the  respective  di¬ 
rections  as  follows: 
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Phase  one  is  then  completed  by  updating  the  internal  energy 
in  the  partial  cell  to  the  intermediate  state  using  the 
same  equation  as  used  for  fixed  boundary  cells  (Equation  28) . 
The  computational  cycle  is  again  completed  by  treat¬ 
ing  mass,  momentum  and  energy  transport  in  phase  two.  In 
this  phase  the  solid  boundary  is  considered  fixed  so  that 
the  treatment  of  the  transport  terms  is  identical  to  that 
described  above  for  the  fixed  boundary  cells. 

For  application  to  the  projectile  boundary,  in  the 
muzzle  blast  calculation,  the  method  described  here  is 
considerably  simplified  since  for  rigid  body,  axisymmetric 
motion; 


0  (t)  =  0  ,  vq(t)  =  vp(t)  (41) 


for  all  beads  on  the  projectile  boundary.  Vp  is  the  pro¬ 
jectile  velocity  whose  time  history  is  determined  by  the 
rigid  body  projectile  dynamics  based  on  the  integrated 
pressure  forces  on  the  projectile  surface.  This  is  des¬ 
cribed  in  Section  4.3  below. 

As  with  the  treatment  of  fixed  boundaries,  the 
numerical  treatment  of  moving  boundaries  given  in  this 
section  is  one  of  numerous  alternative  techniques  which 
could  be  devised.  In  fact,  in  the  present  study,  one 
alternate  technique  was  examined.  In  the  method,  the 
Lagrangian  equation  of  motion  normal  to  the  solid  bound¬ 
ary  was  solved  for  the  pressure  in  the  boundary  cell. 
However,  when  applied  in  the  test  case  *  an  impulsively 
accelerated  projectile  (Section  5.1),  the  method  led  to 
unacceptable  pressure  oscillations  in  the  early  stages 
of  the  calculation.  Thus  it  was  decided  to  use  the  method 
described  above,  in  which  the  pressure  is  not  updated  in 
phase  one  but  is  calculated  at  the  end  of  phase  two  as  in 
a  normal  cell.  The  relatively  good  results  calculated 
for  the  same  test  case  using  the  present  method  lend  some 
justification  to  its  use.  The  test  cases  described  in 
Section  5  indicate  that  the  method  given  here  is  quite 
stable  and  of  acceptable  accuracy,  at  least  for  one  dimen¬ 
sional  rigid  body  boundary  motion. 
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3.4  Stability  and  Accuracy  Considerations 


The  finite  difference  scheme  given  above  is  first 
order  accurate  in  both  mesh  spacing  (Ar,  Az)  and  time  in¬ 
crement  (At).  This  means  that  the  truncation  errors,  in¬ 
troduced  as  a  result  of  the  approximation,  are  linearly 
proportional  to  Ar,  Az  and  At.  Both  accuracy  and  stabil¬ 
ity  of  the  numerical  scheme  are  closely  related,  and  are 
affected  by  the  functional  form  of  the  truncation  errors. 
In  a  time  dependent  calculation,  numerical  instabilities 
can  be  unbounded,  which  are  catastrophic,  or  may  take  the 
form  of  bounded  oscillations  in  flow  quantities.  The  lat¬ 
ter  clearly  affect  the  accuracy  of  the  solution  and  al¬ 
though  they  may  be  within  acceptable  bounds,  they  are 
annoying.  The  overall  accuracy  of  the  time  dependent  so¬ 
lution  depends,  in  addition  to  stability  considerations, 
upon  the  manner  in  which  errors  accumulate  or  cancel  with 
time.  This  is  not  subject  to  analysis,  however,  so  that 
a  final  determination  of  solution  accuracy  must  necessar¬ 
ily  depend  upon  comparisons  to  experiment  or  exact  solu¬ 
tions  . 


Rich1  and  others  have  analyzed  the  nature  of  the 
truncation  errors  implicit  in  the  fluid-in-cell  differed 
cing  technique,  and  have  shown  that  they  take  the  form 
of  viscous  diffusion.  For  example  in  the  axial  momentum 
equation  the  leading  order  truncation  error  is  of  the 
form: 


e 


11  3  /  .  3v\ 

2  F  gJf(Pur4r  -i7) 


+ 


1  _a_ 

2  dz 


(42) 


The  corresponding  terms  in  the  radial  momentum  and  energy 
equations  are  similar.  It  is  noted  that  the  truncation 
error  is  akin  to  a  viscous  diffusion  with  a  non-isotropic 
"numerical  viscosity"  which  depends  on  the  local  flow  ve¬ 
locity,  density  and  mesh  spacing.  The  diffusive  nature 
of  the  truncation  error  terms  is  in  fact  essential  to  the 
success  of  the  fluid-in-cell  technique.  These  terms  ef¬ 
fectively  stabilize  the  difference  equations  to  small  per¬ 
turbations  introduced  by  boundary  conditions.  They  also 
serve  the  important  function  of  smoothing  out  discontinu¬ 
ities  such  as  shock  fronts,  thus  permitting  the  calcula¬ 
tion  of  these  important  phenomena  without  special  treat¬ 
ment. 
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The  effective  numerical  viscosity  is  seen  to  be 
small  in  regions  of  the  flowfield  with  low  density  or 
velocity.  This  is  the  reason  that  oscillations  develop  in 
stagnation  regions  (velocity  0)  or  in  the  rapid  expan¬ 
sion  following  a  strong  shock,  as  noted  by  other  inves¬ 
tigators.  These  oscillations  can  be  reduced  by  the.  in¬ 
troduction  of  an  "artificial  viscosity"  in  the  form  of 
a  viscous  stress.  The  various  forms  possible  are  dis¬ 
cussed  in  standard  works  which  describe  time  dependent 
methods1* .  No  such  artificial  viscosity  is  used  in  the 
present  study,  since  it  is  believed  that  stable  muzzle 
blast  calculations  can  be  performed  without  it.  Thus 
the  additional  inaccuracies  resulting  from  an  artificial 
viscosity  are  avoided. 

The  basic  stability  consideration  in  the  present 
numerical  technique  is  the  Courant  condition.  All  ex¬ 
plicit  difference,  methods  are  subject  to  the  restrictions 
that  the  "domain  of  dependence'"  of  the  finite  difference 
approximation  be  such  as  to  include  the  domain  of  depen¬ 
dence  of  the  physical  process  being  modeled.  In  the 
present  form  it  requires  that  a  signal  cannot  travel  more 
than  one  cell  width  during  a  time  step  At.  For  a  given 
finite  difference  grid  this  places  a  restriction  cn  the 
maximum  allowable  time  step.  In  this  study  the  form  of 
the  Courant  condition  used  is: 


At  < 
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(43) 


where  a  is  the  sound  speed  and  all  cells  in  the  grid  are 
tested.  This  is  the  most  restrictive  form  of  the  Courant 
condition  and  is  valid  for  both  subsonic  and  supersonic 
flow  speeds.  The  restriction  on  At  also  satisfies  the 
requirement  that  an  Euler ian  cell  cannot  be  emptied  dur¬ 
ing  one  time  step.  It  is  noted  in  passing  that  the  Courant 
condition  is  the  basic  factor,  apart  from  the  efficiency 
of  the  overall  scheme,  which  determines  the  computer  time 
required  to  perform  an  unsteady  calculation  to  a  certain 
physical  time.  For  the  muzzle  blast  problem,  the  region 
in  and  around  the  muzzle  is  of  most  interest  and  will 
therefore  be  the  region  in  which  small  cells  are  to  be 
concentrated  for  accurate  resolution.  This  region  of  the 
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flow  is  the  most  restrictive  in  term?  of  the  Courant  con¬ 
dition  since  flow  velocity  and  temperature  (therefore 
sound  speed)  are  greatest  there.  Thus  the  design  of  an 
efficient  grid  for  the  muzzle  blast  problem  will/  as  usual, 
require  a  trade-off  between  solution  accuracy  and  computer 
time. 

As  a  final  remark  it  is  noted  that  the  numerical 
treatment  of  boundary  conditions  can  have  an  important 
effect  on  computational  stability.  This  is  an  area  which 
is  difficult  to  examine  analytically  so  that  it  is  gener¬ 
ally  studied  by  numerical  experimentation.  In  a  qualita¬ 
tive  sense  improper  treatment  of  boundaries  can  generate 
spurious  signals  which  can  be  transmitted  throughout  the 
grid  thereby  affecting  solution  accuracy  and  possibly  sta¬ 
bility,  if  the  signals  become  concentrated  in  low  speed 
regions  of  the  flow.  In  the  muzzle  blast  problem  the 
solid  muzzle  boundaries  and  the  farfield  boundaries  are 
important  in  this  regard. •  The  effect  of  solid  boundaries 
is  discussed  with  regard  to  the  test  case  results  given 
in  Section  5.0.  The  treatment  of  farfield  boundaries  to 
eliminate  spurious  signals  which  could  be  reflected  back 
into  the  flowfield  is  discussed  in  Section  4.0  below.  A 
more  explicit  way  in  which  boundaries  can  affect  stabil¬ 
ity  is  by  way  of  the  reduced  cell  dimensions  of  partial 
boundary  cells.  Partial  cells  which  are  much  smaller 
than  their  full-sized  counterparts  should  thus  be  avoided 
if  possible  when  setting  up  the  grid  and  defining  bead 
locations. 


4.0 


SPECIFICS  OF  METHOD  FOR  MUZZLE  BLAST  PROBLEM 


The  numerical  method  described  above  has  somewhat 
general  application  to  fluid  mechanic  problems  with  fixed 
and  moving  solid  boundaries.  However,  calculations  perfor¬ 
med  to  date  with  the  computer  program,  SAMS,  have  been  rele¬ 
vant  to  the  muzzle  blast  problem.  Some  specifics  of  the 
method  for  application  to  the  muzzle  blast  problem  are 
presented  in  this  section.  The  items  discussed  include: 
muzzle  gas  equation  of  state  in  Section  4.1,  specification 
of  the  time  history  and  transient  boundary  condition  de¬ 
fining  flow  out  of  the  barrel  and  treatment  of  farfield 
boundaries  in  Section  4.2,  and  coupling  of  projectile 
motion  to  the  flowfield  in  Section  4.3. 


4 . 1  Equation  of  State 


The  SAMS  code,  as  presently  configured,  allows 
only  one  equation  of  state  This  is  somewhat  of  a  limi¬ 
tation  for  the  muzzle  blast  problem  since  two  different 
gases  are  involved;  the  high  temperature  and  high  density 
combustion  products  which  make  up  the  muzzle  jet  and  the 
ambient  air.  Fortunately,  both  gases  have  similar  molecu 
lar  weights  and  can  be  expected  to  follow  a  perfect  gas 
law  at  low  fluid  densities.  Since  the  flow  of  muzzle 
gases  in  and  around  the  muzzle  device  and  projectile, 
is  of  most  interest  to  this  study,  it  was  decided  to  use 
an  equation  of  state  characteristic  of  the  muzzle  gas. 

The  effect  which  this  has  on  the  accuracy  of  the  flow  of 
ambient  air  driven  by  the  muzzle  jet  is  commented  on  be¬ 
low. 

The  equation  of  state  used  in  this  study  is  the 
Nobel-Abel  equation  of  state  which  for  our  purposes  can 
be  written  in  the  form: 


(Y-l)  pi 

(1-Bp) 


(44) 


where  y  is  the  ratio  of  specific  heats  and  8  is  the  covol¬ 
ume.  The  values  of  y  and  8  which  are  typical  of  the  muzzle 
gases  and  therefore  used  here  are: 
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Y  =  1.24 


(45) 

0  =  o.ooi  m3/kg 


It  is  noted  that  Equation  (44)  is  the  per fact  gas  law 
with  a  correction  for  Vander  Waal's  forces.  The  equation 
should  be  quite  accurate  for  the  muzzle  gases.  Two  types 
of  errors  are  introduced  in  applying  this  equation  of 
state  to  the  ambient  air.  The  first  is  minor  and  concerns 
the  effect  of  the  covolume  correction.  For  the  strongest 
air  shocks  expected,  p  will  be  considerably  less  than 
10  kg/m3  .  Thus  the  covolume  correction  with  0  =  0.001 
m3/kg  will  introduce  an  error  of  less  than  1.0  percent. 

A  more  appropriate  value  for  the  ratio  of  specific  heats 
for  air  is  y  =  1.4,  so  that  a  more  serious  error  is  intro¬ 
duced  by  using  the  effective  y  of  1.24  for  air.  Thus  the 
pressure  jump  and  velocity  of  the  shock  wave  in  the  ambi¬ 
ent  air  and  the  projectile  blunt  body  shock  will  be  con¬ 
siderably  lower  than  they  should  be. 

Since  the  prime  objective  of  this  study  was  to 
develop  and  verify  the  numerical  technique,  no  attempt 
was  made  to  develop  an  equation  of  state  which  more  accu¬ 
rately  modeled  both  the  muzzle  gases  and  ambient  air. 

A  more  accurate  treatment  could  be  accomplished  in  a  num¬ 
ber  of  ways.  For  example,  the  single  Nobel-Abel  equa¬ 
tion  of  state  could  be  retained  but  with  a  varying  y. 

This  could  be  made  quite  accurate  by  defining  y  in  such 
a  manner  that  it  matches  the  value  for  air  at  low  density 
or  internal  energy  and  matches  values  typical  of  the  muz¬ 
zle  gases  at  high  densities  and  internal  energies.  An 
alternate,  and  potentially  more  accurate,  technique  would 
be  to  use  the  tracer  particles  which  define  the  contact 
surface  between  the  muzzle  gases  and  air,  to  differenti¬ 
ate  between  regions  of  the  flowfield  in  which  the  muzzle 
gas  or  air  equation  of  state  is  to  be  used.  Variations 
of  such  a  technique  are  used  in  later  versions  of  the  OIL 
code  (D0RF10,  HELP11)  to  permit  multi-fluid  calculations. 


4.2  Muzzle  Flow  and  Farfield  Boundary  Conditions 


In  addition  to  the  solid  boundary  condition  dis¬ 
cussed  in  Section  3.3,  two  additional  types  of  boundary 
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conditions  occur  in  the  muzzle  blast  problem.  They  oc¬ 
cur  in  many  flow  problems  and  are  generally  known  as  con¬ 
tinuous  inflow  and  outflow  boundary  conditions,  respec¬ 
tively.  The  inflow  boundary  is  needed  to  define  the  flow 
of  muzzle  gases  into  the  solution  field  at  some  location 
inside  the  gun  barrel,  upstream  of  the  muzzle  device. 

Also  outflow  boundary  conditions  are  defined  on  the  far 
boundaries  of  the  finite  difference  grid.  The  numerical 
treatment  of  both  boundaries  are  discussed  in  this  section. 


Figure  7,  Schematic  of  Numerical  Solution  Domain 


The  muzzle  gas  flow  out  of  the  gun  barrel  is  treated 
as  a  transient  boundary  condition  fixed  at  the  muzzle  exit: 
(j=JB)  as  shown  in  Figure  7  above.  This  boundary  condition 
is  treated  in  a  relatively  straightforward  manner  by  defin¬ 
ing  the  flow  properties  in  fictitious  cells  based  on  curve- 
fits  to  experimental  data  or  independent  calculations  of 
the  flow  out  of  the  muzzle.  Thus,  set 
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The  properties  in  the  fictitious  cells  are  updated  in 
this  manner  in  phase  one  of  the  computational  cycle  and 
used  in  the  appropriate  finite  difference  forms.  They 
are  also  used  in  phase  two  to  determine  the  transport  of 
mass,  momentum  and  energy  through  the  top  boundary  of 
the  cells  at  j  =  JB,  0  <  i  £  IIN,  in  the  usual  donor  cell 
form.  The  time  dependent  properties  (PM,  Pw,  etc. )  for 
an  M-16  rifle  are  given  in  Appendix  B  and  used  in  the 
manner  described  here  for  the  test  case  calculation  pre¬ 
sented  in  Section  5.0. 

An  alternate,  and  more  exact,  method  for  treating 
the  muzzle  flow  properties,  would  be  to  include  the  in¬ 
side  of  gun  barrel  in  the  finite  difference  grid.  Ini¬ 
tial  conditions  (t  =  0)  for  the  gas  inside  the  barrel 
could  be  defined  very  accurately  by  an  internal  ballis¬ 
tics  solution.  The  muzzle  blast  calculation  would  then 
be  initiated  when  the  projectile  moves  to  permit  the  gas 
to  vent.  This  would  result  in  an  expansion  wave  which 
would  proceed  upstream  through  the  column  of  gas  in  the 
barrel.  The  finite  difference  grid  inside  the  barrel 
could  be  quite  coarse  in  comparison  to  the  grid  just  out¬ 
side  the  barrel.  In  this  method,  the  flow  inside  the 
barrel  would  be  coupled  to  the  resulting  muzzle  blast 
flowfield  so  that  no  approximation  (other  than  finite 
difference)  would  be  involved.  On  the  time  scale  of 
transitional  ballistics  (t  <  100  microseconds)  the  muzzle 
flow  properties  vary  but  a  few  percent  so  that  either  of 
the  above  methods  should  be  quite  accurate  for  the  pur¬ 
pose  of  this  study. 

The  other  boundary  condition  needed  in  the  muzzle 
blast  computation,  involves  the  far  boundaries  of  the 
finite  difference  grid  as  shown  in  Figure  7.  The 
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axis  of  symmetry  (r  =  0)  is  treated  as  a  reflective  bound¬ 
ary  as  already  discussed.  The  other  boundaries  of  the 
grid  (j  =  0,  j  =  JMAX  and  i  =  IMAX)  must  also  be  treated 
in  some  sense.  It  is  noted  that  the  interaction  of  the 
flow  with  grid  boundaries  is  almost  always  a  problem  in 
finite  difference  calculations.  If  the  boundary  is  not 
treated  in  a  physical  manner,  spurious  signals  can  be 
generated  which  may  destroy  a  time  dependent  calculation. 
The  most  difficult  type  of  interaction  to  treat  is  the 
passage  of  a  shock  through  the  boundary  which  unfortu¬ 
nately  could  occur  in  the  muzzle  blast  problem.  Ideally 
the  boundary  would  completely  absorb  (or  transmit)  the 
shock  but  in  practice  this  is  difficult  to  accomplish 
since  the  flow  outside  the  grid  clearly  depends  in  a  com¬ 
plex  manner  on  the  flow  inside  the  grid.  With  this  prob¬ 
lem  ir.  mind,  it  is  a  good  idea  to  place  the  far  boundaries 
far  enough  away  from  the  flow  region  of  interest  (muzzle) 
so  that  the  calculation  either  terminates  before  the  shock 
reaches  the  boundary  or  before  any  reflected  signals  can 
interact  with  the  flow  regions  of  interest.  This  can  be 
accomplished  by  stretching  the  grid. 

In  any  event  the  method  presently  used  for  these 
boundaries  is  known  as  the  "transmittive"  or  continuous 
"outflow"  boundary  condition4.  The  method  simply  involves 
defining  flow  properties  in  fictitious  cells  outside  tne 
boundary  to  be  equal  to  the  flow  properties  in  the  adja¬ 
cent  active  cells  inside  the  boundary.  This  is  done  at 
all  stages  of  the  calculation.  This  treatment  ir  essence 
fixes  the  flow  gradients  equal  to  zero  at  the  grid  bound¬ 
aries,  and  clearly  can  be  accurate  only  very  far  from  any 
active  region  of  the  flowfield. 


4.3  Dynamics  of  Projectile  Motion 


As  noted  in  an  earlier  section,  the  projectile  is 
initially  accelerating  with  respect  to  the  muzzle  due  to 
the  high  muzzle  gas  pressure  at  its  base.  At  later  times 
it  decelerates  as  the  muzzle  gas  expands  around  the  pro¬ 
jectile  thereby  decreasing  the  base  pressure.  In  this 
study,  the  response  of  the  projectile  to  the  time  vary¬ 
ing  pressure  forces,  exerted  on  its  surface, is  treated. 
Since  the  calculation  is  inviscid,  only  pressure  forces 
are  relevant.  Also  the  calculation  is  axisymmetric,  so 
that  only  the  axial  motion  of  the  projectile  need  be  con¬ 
sidered. 


\ 
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The  projectile  equation  of  motion  in  the  axial 
direction  is  simply: 


dvn 

Mp  ~dt 


■jC 


p(g 


ft)dA 


(47) 


where  Mp/  Vp  are  the  mass  and  velocity  of  the  projectile 
respectively  and  £z  is  the  unit  vector  in  the  axial  di¬ 
rection  and  ft  is  the  unit  normal  to  the  projectile  sur¬ 
face  element  dA.  The  finite  difference  approximation  to 
Equation  (47)  is  used  to  update  the  projectile  velocity 
to  the  advanced  time  tn+^  starting  with  its  initial  muz¬ 
zle  velocity  at  time  t  =  0.  Summing  the  pressure  forces 
in  each  of  the  partial  cells  on  the  projectile  boundary, 
results  in  the  following  finite  difference  form: 
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where  the  summation  is  over  all  the  boundary  cells  with 

index  k  . 

P 

Over  a  time  scale  typical  of  transitional  ballis¬ 
tics  (t  <  100  microseconds)  the  change  in  projectile  velo' 
city  is  negligible.  In  calculations  performed  to  date, 
the  change  has  been  less  than  1  m/sec.  As  a  result,  the 
SAMS  code  has  a  constant  projectile  velocity  option  which 
eliminates  the  need  for  performing  the  calculation  des¬ 
cribed  here, 
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5.0  TEST  CASE  RESULTS 


As  has  been  stated  a  number  of  times  in  this  re¬ 
port,  the  main  innovation  in  the  numerical  technique  pre¬ 
sented  here  is  the  treatment  of  fixed  or  moving  solid 
boundaries  in  an  Eulerian  grid.  It  is  well  known  that 
boundary  conditions  are  often  crucial  to  any  finite  dif¬ 
ference  technique  for  fluid  mechanics  problems.  This  is 
especially  true  of  the  muzzle  blast  problem  since  so  many 
of  the  interesting  phenomena  are  a  result  of  flow  inter¬ 
action  with  fixed  or  moving  boundaries.  Thus  the  first- 
test  case,  presented  in  Section  5.1,  was  performed  to 
verify  the  accuracy  and  stability  of  the  treatment  of 
moving  boundaries.  The  complete  computer  program  capa¬ 
bility  is  exercised  in  the  typical  muzzle  blast  test  cases 
presented  in  Sections  5.2  and  5.3.  These  calculations  sim¬ 
ulate  the  muzzle  gas  flow  of  an  M-16  rifle  with  and  *?ith- 
out  a  simple  muzzle  device 


5 . 1  Blunt  Projectile  Impulsively  Started  From  Rest 


This  test  case  involves  the  calculation  of  the  flow 
around  a  flat-faced  cylindrical  projectile  accelerated  im¬ 
pulsively  to  1500  m/sec  in  an  ideal  gas  (y  =  1.3)  at  rest. 
The  calculation  was  performed  by  passing  the  rigid  bound¬ 
ary  defining  the  projectile  surface  through  a  fixed  Eu¬ 
lerian  grid.  The  purpose  of  the  calculation  was  to  check 
out  the  logic,  accuracy  and  stability  of  the  partial  cell 
technique  described  above  for  moving  solid  boundaries. 

A  schematic,  describing  the  geometry  and  parameters 
for  this  test  case,  is  presented  in  Figure  8  below.  The 
cylinder  has  a  radius  of  0.5  cm  and  length  of  0.5  cm  and 
its  velocity  is  taken  as  1500  m/sec.  The  ambient  gas  is 
considered  ideal  with  y  =  1.3.  For  these  conditions  the 
projectile  Mach  number  is  M*,  =  5.0.  The  finite  differ¬ 
ence  grid  was  made  up  of  zones  with  dimensions  Ax  =  Ay  = 
0.1  cm,  and  its  extent  was  1.5  cm  in  radius  by  2.0  cm  in 
the  axial  direction.  The  maximum  allowable  time  step  was 
determined  by  the  Courant  condition  with  Courant  number  = 
0.25.  No  hint  of  numerical  instability  was  evident  in 
the  calculated  results,  indicating  that  a  larger  Courant 
number  could  have  been  used  with  success. 
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Theoretically,  the  axis  of  the  cylinder  should  reach  these 
conditions  at  steady  state.  Figures  9,  10,  11  given  below, 
give  the  time  history  of  the  pressure,  temperature  and 
density  respectively,  of  the  finite  difference  zone  on 
the  axis,  just  ahead  of  the  projectile.  These  results  are 
shown  compared  to  the  theoretical  stagnation  conditions 
and  the  comparison  is  judged  good  considering  the  crude¬ 
ness  of  the  mesh.  The  over-prediction  of  temperature  is 
consistent  with  the  dissipative  nature  of  the  higher  order 
terms  resulting  from  the  discretization  process. 

Experimental  results  for  shock  stand-off  distance13 
for  a  flat-faced  cylinder  at  =  5  indicate  that: 


(49) 


Thus  for  the  1  cm  diameter  projectile  of  this  test  case 
As  =  0.3  cm.  The  numerical  results  are  consistent  with 
this  experimental  result.  If  the  shock  is  defined  by  the 
maximum  pressure  or  density  gradient,  the  numerical  re¬ 
sults  predict  a  shock  in  the  third  finite  difference  zone 
in  front  of  the  projectile  or  As| numerical  =  0*25.  This 
compares  to  the  numerical  result  within  the  accuracy  of 
the  finite  difference  mesh. 

One  final  check  on  the  computed  results  can  be  ob¬ 
tained  by  comparing  the  radial  velocity  on  the  face  of  the 
cylinder  with  theoretical  predictions.  The  constant  den¬ 
sity  shock  layer  solution  of  Probstein1  ''  predicts  a  radial 
velocity  distribution  given  by: 


where  Uro  or  Vp  is  the  velocity  of  projectile,  ps  and  pro 
are  the  stagnation  density  and  ambient  density,  respectively. 


As  is  the  shock  stand-off  distance  and  r,  R^  are  the  radius 
and  body  radius,  respectively.  This  solution  compares 
favorably  with  experiment  for  r/Rfc  a  0.5.  Comparison  of 
the  numerical  results  (at  t  -  7  microseconds)  with  this 
theoretical  radial  velocity  distribution  are  shown  in 
Figure  12  below.  The  comparison  is  again  judged  to  be 
quite  good.  The  fact  that  the  numerical  result  is  below 
the  theoretical  indicates  that  the  calculation  may  not 
have  reached  steady  state. 

This  test  case  is  a  relatively  severe  test  of  the 
numerical  treatment  of  moving  boundaries.  The  comparison 
of  results  with  experiment  and  theory  is  good  and  within 
the  accuracy  of  the  relatively  crude  finite  difference 
mesh.  The  mesh  was  chosen  as  typical  of  the  mesh  size 
used  in  the  region  of  the  projectile  in  a  complete  muzzle 
blast  calculation  considering  practical  usage  of  computer 
storage  and  computational  time.  No  hint  of  numerical  in¬ 
stability  or  bounded  oscillations  was  evident  in  the  cal¬ 
culation.  The  only  anomaly  was  the  nonuniformity  of  the 
time  history  cf  flow  pressure  and  density  in  the  partial 
boundary  cells  as  shown  in  Figures  9  and  11.  This  is  not 
serious,  however,  and  it  is  believed  that  finer  zoning 
would,  result  in  smoother  results. 
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Figure  13.  Schematic  of  M-16  Muzzle  Blast  Calculation 


The  SAMS  code  was  used  to  calculate  the  muzzle 
blast  flowfield  of  an  M-16  rifle  and  the  calculated  re¬ 
sults  are  presented  in  this  section.  The  muzzle  and  pro¬ 
jectile  geometry  and  finite  difference  grid  design  are 


shown  schematically  in  Figure  13.  The  muzzle  consists 
of  a  cylinder  with  an  inside  diameter  of  0.556  cm  (5.56 
mm  bore)  and  outside  diameter  of  1.556  cm.  The  Nobel- 
Abel  equation  of  state ,  described  earlier,  was  used  for 
both  the  muzzle  gas  and  ambient  air  and  the  transient  flow 
properties  at  the  muzzle  exit  were  determined  from  an  in¬ 
dependent  calculation  of  Celmins*2  (see  Appendix  B) . 

Based  on  these  flow  properties,  a  projectile  muzzle  velo¬ 
city  of  1210  m/sec  was  used  and  was  assumed  constant 
throughout  the  calculation.  The  finite  difference  grid 
was  varied  in  a  continuous  manner  from  Ar  =  0.055  cm, 

Az  =  0.075  cm  near  the  muzzle  exit  to  Ar  =  Az  =  0.5  cm 
near  the  boundaries  of  the  grid.  The  projectile  was  de¬ 
fined  by  a  sequence  of  "beads,"  with  its  base  initially 
(t  =  0)  at  the  muzzle  exit.  Finally  the  integration  time 
step  was  determined  from  the  Courant  condition  with  a 
Courant  number  of  0.3. 

Computer  generated  plots  of  the  calculated  results 
are  presented  in  Figures  14  through  27  at  the  end  of  this 
section.  Figures  14  through  18  are  velocity  vector  plots 
of  the  flowfield  at  increments  of  about  5  microseconds 
(y  sec)  (every  100  cycles)  from  t  =  5.91  y  sec  (cycle  100) 
to  t  =  25.18  y  sec  (cycle  500).  Figures  19  through  23  are 
pressure  contour  plots  for  the  same  times  as  the  velocity 
vector  plots.  Pressure  contours  of  2,  10  and  100  bars  are 
plotted.  Finally,  Figures  24  and  25  are  density  contour  and 
temperature  contour  plots  respectively  at  t  =  10.93  y  sec 
(cycle  200)  and  Figures  26  and  27  are  corresponding  plots  at 
t  =  25.18  y  sec  (cycle  500).  Contours  of  0.5,  2.0,  5.0  and 
10.0  kg/ra3  are  plotted  in  the  density  plots  and  contours 
of  500,  1,000,  2,000  and  3,000°K  are  given  in  the  temper¬ 
ature  plots. 

The  figures  clearly  show  the  development  of  all 
the  expected  phenomena.  The  advancing  shock  wave  in  the 
ambient  air  which  weakens  to  a  compression  wave  as  it 
diffracts  around  the  gun  barrel  is  shown  most  clearly  in 
the  pressure  contour  plots  (Figures  19  through  23) <  The 
transient  development  of  the  centered  expansion  wave  origi¬ 
nating  at  the  corner  of  the  gun  barrel  is  also  shown  in 
the  pressure  plots.  The  relatively  sharp  pressure  rise 
through  the  shock  and  its  decay  to  a  level  below  ambient 
with  subsequent  rise  in  the  muzzle  jet  can  also  be  seen. 

This  behavior  is  characteristic  of  a  spherical  blast  wave 
and  is  predicted  by  the  calculation.  Also  note  the  devel¬ 
opment  of  the  projectile  blunt  body  shock  and  its  interac¬ 
tion  with  the  air  shock. 
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The  development  of  the  jet  of  muzzle  gas,  shown 
in  both  velocity  vector  and  pressure  plots  is  seen  to 
be  characteristic  of  an  under-expanded  jet.  It  is  noted 
that  the  jet  soon  reaches  a  "quasi-steady"  state  in  the 
region  around  the  exit.  This  is  evidenced  by  the  stable 
position  of  the  10  bar  pressure  contour  from  5.91  p  sec 
on,  and  the  stable  position  of  the  2  bar  contour  from  20 
jj  sec  on.  The  relatively  steady  muzzle  jet  flow  is  also 
shown  in  the  density  and  temperature  contour  plots  (Fig¬ 
ures  24  through  27) . 

Two  phenomena  shown  in  the  results  are  of  partic¬ 
ular  note.  The  first  is  the  separation  of  the  flow  at  the 
corner  of  the  muzzle  face  as  shown  in  the  velocity  vector 
plots  from  10.93  p  sec  to  25.18  p  sec  (Figures  15  through  18). 
The  resulting  recirculating  flow  region  and  mixing  layer 
which  extends  from  the  corner  and  perpendicular  to  the 
barrel  are  clearly  indicated.  The  flow  seperation  is  a 
physically  real  phenomena,  not  to  be  expected  in  an  invis- 
cid  calculation,  and  is  due  to  the  numerical  viscosity  in 
the  finite  difference  scheme.  The  second  result  of  note 
is  the  shock  wave  in  the  muzzle  jet  at  the  base,  of  the  pro¬ 
jectile.  It  is  most  clearly  shown  in  the  pressure,  density 
and  temperature  plots  at  t  =  10.9  p  sec  (Figures  20,  24 
and  25  respectively) .  This  is  due  to  the  expansion  of 
the  muzzle  jet  to  a  supersonic  velocity  relative  to  the 
projectile,  and  provides  additional  evidence  that  the  nu¬ 
merical  treatment  of  moving  solid  boundaries,  developed 
in  this  study,  is  accurate. 

The  results  presented  here  are  physically  reasonable 
in  every  way  and  are  indicative  of  the  quality  of  results 
which  can  be  obtained  with  the  SAMS  code.  Although  no  er¬ 
ror  estimates  or  extensive  experimental  comparisons  have 
been  performed,  the  results  are  believed  to  be  accurate. 

All  flow  quantities  are  smoothly  varying  from  cell  to  cell 
and  no  hint  of  numerical  oscillation  or  instability  is  in¬ 
dicated  with  one  exception.  The  exception  is  the  existence 
of  bounded  oscillations  in  flow  quantities  in  an  isolated 
region  of  the  flow  on  the  low  velocity  side  of  the  mixing 
layer  described  above.  The  oscillations  are  evident  in 
the  pressure  plots,  starting  at  cycle  300  (Figure  21)  and 
extending  to  the  end  of  the  calculation.  They  do  not  des¬ 
troy  the  calculation  and  remain  in  that  region  of  the  flow. 

It  is  noted  that  they  occur  in  a  low  velocity  (<  200  m/sec) 
and  low  density  (<  0.5  kg/m3)  region  of  the  flow.  It  is  re¬ 
called  that  the  numerical  viscosity  in  the  finite  difference 
scheme  becomes  ineffective  for  low  velocities  and  low  densi¬ 
ties  so  that  the  oscillations  could  be  due  to  the  lack  of 
numerical  smoothing.  It  is  important  to  point  out  that  the 


oscillations  occur  in  the  low  velocity  region  of  the  mixing 
layer.  It  is  well  known  that  a  mixing  layer  separating  a 
low  velocity  from  a  high  velocity  flow  is  unconditionally 
unstable15/  developing  in  turn  vortex  oscillations  and  ul¬ 
timately  turbulence.  Thus  the  oscillations  could  be  a  re¬ 
sult  of  the  attempt  of  the  inviscid  theory  to  model  a  real 
physical  phenomena „  It  is  believed  that  the  introduction 
of  an  artificial  viscosity  into  the  numerical  technique 
wculd  damp  out  the  oscillations  but  this  was  not  attempted 
in  this  study.  An  artificial  viscosity  could  have  adverse 
effects  on  other  regions  of  the  flow  so  that  its  use  should 
be  accompanied  by  extensive  numerical  experimentation. 

One  final  comment  which  is  of  interest  is  that  the 
calculation  described  in  this  section  was  performed  in  less 
than  16  minutes  of  CPU  time  on  a  CDC  6600.  This  indicates 
that  muzzle  blast  calculations  can  be  performed  with  com¬ 
puter  program  SAMS  with  a  relatively  modest  investment  in 
computer  time . 
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Figure  23,  Pressure  Contour  Plot:  M-16  Rifle 
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Figure  28.  Schematic  of  M-16  with  Muzzle  Device 


In  this  section,  the  results  of  a  SAMS  calculation 
of  the  muzzle  blast  flowfield  of  an  M-16  rifle  with  a  sim¬ 
ple  muzzle  device  are  presented.  The  muzzle,  muzzle  device 
and  projectile  geometry  and  finite  difference  grid  design 
are  shown  schematically  in  Figure  28.  The  barrel  and 


projectile  are  identical  to  those  used  in  the  calcula¬ 
tion  of  the  previous  section.  The  muz2le  device  consists 
of  an  extension  of  the  barrel  with  circumferential  slits 
to  permit  the  flow  to  expand  radially.  The  equation  of 
state,  muzzle  jet  time  history,  muzzle  velocity  and  basic 
finite  difference  grid,  used  in  this  calculation,  are 
identical  to  those  used  in  the  M-16  calculation  presented 
above.  The  two  closed  solid  boundaries  of  the  muzzle  "baf¬ 
fles"  were  defined  by  a  sequence  of  fixed  beads  introduced 
into  the  solution  field  at  the  positions  shown  in  the  fig¬ 
ure.  The  resulting  fixed  partial  cells  along  the  solid 
boundaries  were  just  slightly  smaller  than  the  correspond¬ 
ing  rectangular  cells  in  the  region. 

Computer  generated  plots  of  the  calculated  results 
are  shown  in  Figures  29  through  42,  at  the  end  of  this  sec¬ 
tion.  Figures  29  through  33  are  velocity  vector  plots  of 
the  flowfield  at  increments  of  about  5  microseconds  (y  sec) 
(every  100  cycles)  from  t  =  5.91  y  sec  (cyle  100)  to  t  = 
25.77  y  sec  (cycle  500).  Figures  34  through  38  are  pres¬ 
sure  contour  plots  for  corresponding  times.  Again  pressure 
contours  of  2,  10  and  100  bars  are  plotted.  Finally,  Fig¬ 
ures  39  and  40  are  density  and  temperature  contour  plots 
at  t  =  10.95  y  sec  and  Figures  41  and  42  are  the  corres¬ 
ponding  plots  at  t  =  25.77  y  sec. 

The  figures  show  the  details  of  the  time  develop¬ 
ment  of  the  flow  which  in  this  case  is  much  more  compli¬ 
cated  than  the  basic  M-16  calculation  of  the  previous 
section.  The  most  notable  feature  of  the  flow  is  the  ex¬ 
pected  strong  radial  jet  of  flow  between  the  barrel  and 
the  first  baffle.  The  flow  is  turned  by  a  strong  shock 
resting  on  the  inside  corner  of  the  baffle  and  which  ex¬ 
tends  into  the  jet  and  ultimately  weakens  to  a  compression, 
due  to  its  interaction  with  expansion  waves  from  the  corner 
of  the  barrel  and  the  outside  corner  of  the  baffle.  As  in 
the  previous  calculation,  the  flow  separates  at  the  face  of 
tile  barrel  and  also  at  the  outside  corner  of  the  baffle. 
Based  on  the  strong  radial  diversion  of  the  flow  in  this 
calculation  as  compared  to  the  previous  calculation,  it  is 
expected  that  the  simple  muzzle  device  used  here  would  sig¬ 
nificantly  decrease  recoil  forces  on  the  barrel. 

The  velocity  vector  plots  also  show  the  development 
of  a  second  radial  jet  flowing  through  the  space  between 
the  first  and  second  baffle.  The  jet  begins  to  develop  at 
about  15  y  sec  and  is  considerably  weaker  than  the  jet  des¬ 
cribed  above.  Also  of  note  is  the  flow  inside  the  muzzle 
device  which  is  seen  to  be  very  close  to  one  dimensional 
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throughout  the  calculation.  The  jet  of  muzzle  gases  ex¬ 
pands  to  about  10  bars  then  recompresses  behind  the  pro¬ 
jectile.  The  flow  in  the  region  around  the  barrel  exit 
and  inside  the  first  baffle  is  shown  to  be  approximately 
steady  from  about  15  y  sec  to  the  end  of  the  calculation. 
Also  of  note  is  the  fact  that  the  pressure  at  the  base  of 
the  projectile  remains  above  10  bars  throughout  the  calcu¬ 
lation,  whereas  in  the  previous  calculation  the  base  pres¬ 
sure  was  less  than  10  bars  from  15  y  sec  on.  Thus  the 
projectile  is  continuing  to  accelerate  through  the  muzzle 
device. 


As  with  the  previous  results,  the  results  of  this 
calculation  are  physically  reasonable  and  numerically 
smooth  with  the  exception  of  the  localized  flow  oscilla¬ 
tions  which  develop  in  the  same  region  (mixing  layer  at 
barrel  exit)  as  described  in  Section  5.2.  A  particularly 
pleasing  result  is  the  flow  around  the  baffles  which  again 
is  physically  reasonable  and  numerically  smooth.  This 
provides  some  evidence  that  the  treatment  of  fixed  solid 
boundaries,  developed  in  this  study,  is  accurate.  Finally, 
it  is  noted  that  the  present  calculation  was  also  performed 
in  less  than  16  minutes  of  CPU  time  on  a  CDC  6600. 
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CONCLUSIONS  AND  RECOMMENDATIONS 


The  main  result  of  this  study  is  the  generaliza¬ 
tion  of  the  fluid-in-cell  finite  difference  technique  to 
account  for  arbitrary  fixed  and  moving  boundaries  in  the 
Eulerian  grid.  The  technique  has  been  developed  in  gen¬ 
eral  but  applied  specifically  to  the  muzzle  blast  flow- 
field  problem  in  the  SAMS  computer  program.  The  results 
presented  in  this  report  indicate  that  the  treatment  of 
fixed  and  moving  boundaries  is  valid  and  that  the  overall 
numerical  technique  is  a  good  approach  to  the  simulation 
of  muzzle  blast  flowfields  which  include  a  moving  projec¬ 
tile  and  fixed  muzzle  device  boundary.  The  results  indicate 
that  the  SAMS  code  could  be  a  useful  tool  for  muzzle  de¬ 
vice  studies  since  it  is  relatively  easy  to  set  up  and 
run  and  it  makes  relatively  efficient  use  of  computer 
time.  Useful  muzzle  flow  simulations  can  be  performed 
in  approximately  15  minutes  of  CPU  time  on  a  CDC  6600. 

The  numerical  technique  developed  during  the  study 
could  be  improved  in  a  number  of  ways,  some  of  which  have 
been  mentioned  in  previous  sections.  For  example,  the 
treatment  of  boundary  conditions  could  be  improved.  The 
displaced  center  of  mass  of  the  partial  boundary  cells 
could  be  accounted  for  in  the  Lagrangian  phase  o'f  the 
computational  cycle  and  also  in  extrapolating  flow  quant¬ 
ities  to  the  solid  boundary.  Also  an  improved  treatment 
of  the  farfield  boundaries  to  absorb  the  air  shock  and 
trailing  flow  with  minimum  reflection  would  allow  moving 
the  computational  boundaries  closer  to  the  muzzle.  This 
would  result  in  more  accurate  resolution  of  the  flow  in 
the  muzzle  device,  through  an  increased  concentration  of 
grid  points  near  the  muzzle.  In  a  related  vein,  the  accu¬ 
racy,  stability  and  efficiency  of  the  overall  numerical 
technique  should  be  investigated  to  a  greater  degree  then 
was  possible  in  this  study.  Accuracy  could  be  verified 
by  comparison  to  experimental  data  or  exact  solutions  for 
special  cases  (1-D  projectile  motion,  etc.).  The  stabil¬ 
ity  of  the  technique  could  be  enhanced  by  the  use  of  an 
artificial  viscosity  to  damp  out  flow  oscillations  such 
as  those  which  develop  in  the  unstable  mixing  layer  dis¬ 
cussed  in  the  test  case  calculations.  Finally,  the  effi¬ 
ciency  of  the  numerical  technique  for  application  to  muzzle 
blast  problems  could  be  "fine  tuned"  through  more  exten¬ 
sive  experimentation  *ith  grid  spacing  and  distribution 
and  Courant  numbers.  All  of  the  above  developments  would 
require  numerical  experimentation  to  optimize  and  build 
confidence  in  the  numerical  technique. 
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Certain  generalizations  to  the  present  technique 
would  both  enhance  the  simulation  and  permit  the  study 
of  more  detailed  phenomena.  As  noted  in  Section  4.1, 
the  present  single  equation  of  state  can  be  made  quite 
accurate  for  describing  both  the  muzzle  gas  and  air. 
However,  the  addition  of  a  two  equation  of  state  capa¬ 
bility  would  increase  the  accuracy  and  also  permit  a 
more  detailed  study  of  the  mixing  of  the  muzzle  gases 
and  air.  For  modern  high  velocity  weapons,  two  phase 
flow  effects  in  the  muzzle  jet  can  be  quite  important. 

That  is,  burning  propellent  particles  will  very  likely 
accompany  the  combustion  products  out  of  the  muzzle. 

Such  effects  could  be  studied  by  introducing  particles 
into  the  flow  which  couple  to  the  flow  through  mass  (pro¬ 
ducts)  ,  momentum  (drag)  and  energy  (heat  deposition)  ex¬ 
change.  This  could  be  accomplished  to  varying  degrees 
of  accuracy,  as  is  commonly  done  in  two  phase  flow  compu¬ 
tations  . 

In  summary,  it  is  believed  that  the  present  numer¬ 
ical  technique,  incorporated  in  the  SAMS  code,  can  be  quite 
useful  in  the  study  of  muzzle  gas  flows.  Certain  improve¬ 
ments  to  and  generalizations  of  the  technique,  discussed 
here,  would  enhance  both  its  utility,  accuracy  and  flexi¬ 
bility. 
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APPENDIX  A 


GEOMETRICAL  DESCRIPTION  OF  PARTIAL  CELLS 


The  numerical  technique  for  treatment  of  solid 
boundaries  within  the  Euler ian  grid,  given  above,  re¬ 
quires  certain  geometrical  factors  which  define  the  par¬ 
tial  cells  formed  by  intersections  of  the  boundary  curve 
with  the  rectangular  grid  lines.  These  factors  are  the 
fractions  of  the  cell  sides  open  to  the  flow  (Ak,  k  =  1, 
2,  3,  4) ,  the  fraction-  of  the  rectangular  cell  volume  (F) 
occupied  by  fluid  and  an  angle  (0')  defining  the  local 
normal  to  the  boundary.  The  details  of  the  procedure 
for  calculating  these  factors  are  given  in  this  appendix. 
The  discussion  applies  for  any  fixed  boundary,  or  for  a 
moving  boundary  defined  at  any  time  (t)  by  the  locations 
of  a  sequence  of  boundary  points. 


Figure  A.l.  Schematic  of  Partial  Cell 


A  general  fixed  or  moving  boundary  is  described 
by  a  string  of  beads  fixed  to  the  boundary  and  numbered 
from  1  to  N  in  ascending  order  in  a  counterclockwise 
fashion.  Only  one  bead  is  allowed  to  reside  in  a  fixed 
Eulerian  cell.  In  the  event  of  more  than  one  bead  lo¬ 
cated  in  a  cell,  the  leading  bead  is  used  to  represent 
the  boundary.  Hence  the  density  of  beads  is  related  di' 
rectly  to  the  spatial  resolution  of  the  numerical  solu¬ 
tion. 


Each  cell  i r  defined  by  the  position  of  the  four 
corners  as  indicated  in  Figure  A.l.  In  other  words,  the 
coordinates  of  (ri,Zj),  (ri_i,Zj),  ,  (r^Zj^) 

are  known  for  each  cell.  The  four  sides  of  each  cell  are 
labeled  1,  2,  3  and  4  in  a  counterclockwise  manner  as  al¬ 
so  depicted  in  the  figure.  The  partial  cell  occupied  by 
bead  q  is  characterized  by  five  quantities:  A^(q),  A2(q) 
A3(q),  A4 (q)  which  signify  the  percentage  of  open  area 
at  surface  1,  2,  3  and  4,  respectively,  and  F(q)  repre¬ 
senting  the  percentage  of  the  volume  open.  For  the  pur¬ 
pose  of  this  appendix,  these  quantities  are  associated 
with  bead  q  which  resides  in  the  partial  cell  being  con¬ 
sidered. 

In  order  to  arrive  at  the  above  set  of  quantities, 
we  must  have  some  additional  information  for  each  partial 
cell.  The  nature  of  the  partial  cell  is  characterized  by 
the  two  indices  that  denote  the  two  partially  open  sur¬ 
faces.  For  example,  the  partial  cell  occupied  by  bead  q 
in  Figure  A.l  is  type  1-2  cell,  because  surfaces  1  and  2 
are  partially  opened  but  surfaces  3  and  4  are  completely 
open.  Hence  A^q)  <  1,  A2{q)  <  1  and  A3  (q)  =  A4(q)  =  1 
while  the  percentage  of  the  volume  open,  F{q),  is  less 
than  unity.  One  can  readily  count  that  there  are  six 
types  of  partial  cells:  1-2,  1-3,  1-4,  2-3,  2-4,  and 
3-4. 


Figure  A, 2*  Indexing  of  Partial  Cell  Boundary  Points 
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To  evaluate  A].{q),  ...  A4  (q) ,  F{q),  we  must  have 
the  coordinates  of  the  five  points  marked  1,  2,  3,  4,  5 
shown  in  Figure  A. 2  above.  Namely#  we  must  have  r! (q) , 

(q)  ,  ...  r5  (q) ,  zs(q)  .  Point  1  always  refers  to 
bead  q,  point  2  signifies  the  intersection  of  the  bound¬ 
ary  with  the  surface  denoted  by  the  leading  index,  point 
5  represents  the  intersection  of  the  boundary  with  the 
surface  denoted  by  the  trailing  index.  It  should  be  noted 
that  points  3  and  4  might  coincide,  depending  on  the 
character  of  the  partial  cell,  for  example,  in  the  case 
of  the  type  1-2  cell. 


The  procedure  for  classifying  the  partial  cells 
at  any  time,  t,  depends  on  the  relative  magnitudes  of 
i(q+l),  j (q+1) ,  i(q)  and  j (q) .  For  example,  take  the 
following  case:  i(q+l)  >  i(q),  j (q+1)  =  j (q)  as  shown 
in  Figure  A. 3. 


The  line  joining  beads  q  and  q+1  is  of  the  form: 


z  =  Fz(r;  z(q),  r(q),  z(q+l),  r(q+l))  (A.l) 
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where  the  function  F  is  given  by: 

z 


F 

z 


2  (q)-z  (q+1) 
r  (q)-r(q+l) 


r  + 


[z(q) 


2  (q)-z  (q+1) 
r (q)-r (q+1) 


r  (q)] 


(A. 2) 


then 


z4  =  z[j(q)J  -  2d 


(A.  3) 


where 


ZD  "  P2  (<3)  I  ;  2  (q)  ,  r(q),  z(q+l),  r  (q+1) ) 


and  the  fractional  areas  are: 


A4  (q)  =1 - —  =  A-  (q+1) 

Az(q)  Z 


(A.  4) 


with  Az(q)  signifying  the  size  of  the  cell  in  the  z- 
direction.  Coordinates  of  the  intersection  with  the 
cell  side  are: 


r2(q) 


r5 (q+1) 


z2(q)  =  zD  =  z5 (q+1) 


(A.  5) 
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Bead  q+1,  and  all  successive  beads,  moving  counter¬ 
clockwise  around  the  boundary,  are  treated  in  the  same  man¬ 
ner.  To  sum  up  the  computational  procedure  in  graphical 
form,  without  the  detailed  formulae.  Table  A.l  is  presented 
at  the  end  of  this  section. 

After  the  preceding  treatment  has  been  applied  to 
each  bead  on  the  boundary,  the  leading  (first)  index  for 
bead  q,  as  well  as  the  trailing  (second)  index  for  the 
bead  q+1  is  known.  The  location  of  the  leading  cut  for 
bead  q  and  the  location  of  the  trailing  cut  for  bead  q+1 
are  also  known.  Each  partial  cell  is  assigned  two  indices, 
two  of  the  four  A^(q)'s,  k  =  1,  2,  3,  4,  and  the  location 
of  points  1,  2  and  5.  In  order  to  complete  the  descrip¬ 
tion  of  the  partial  cell,  the  remaining  two  Afc(q) '  s,  k  = 

1,  2,  3,  4  must  be  obtained,  and  the  fraction  of  the 
volume  open,  F(q),  must  be  computed.  That  is,  the  loca¬ 
tions  of  points  3  and  4  for  each  cell  must  be  given.  This 
is  done  according  to  Table  A. 2  at  the  end  of  this  section. 

Once  the  position  of  the  five  key  points  are  found, 
namely,  r  ^(q) ,  z  (q)  ...  rg(q),Ze(q)  ,  the  percentage 

of  the  volume  open  F(q)  can  Be  obtained  from  the  following 
formula . 


F(q) 


(1  -  f  (q)  if  f  (q)  >  0 
U  f(q)  if  f (q)  <  0 


where 


(A.  6) 


,  A  Lzk+i  (<*)“2k(q)]  Lrk+i  +rk  (q)J 2 

f(q)  =  i  - — - 

Lr  (q)  >  ~r‘  (i  (q)  -D]  [2  ( j  (q) )  -z  ( j  (q)  -1)] 


The  expression  f(q)  can  easily  be  obtained  by  integrating 
the  cylindrical  volume  enclosed  by  the  five  points.  To 
complete  the  geometrical  definition  cf  the  partial  bound¬ 
ary  cells,  the  direction  cosines  of  the  outward  directed 
normal  to  the  solid  boundary  is  defined  in  the  following 
approximate  manner  (see  Figure  A. 4). 
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(sin  0')  = 


z2(q)-z5(q) 


>2  (q'  ~z5  (q)  3  +Cr2  (q)  "r5 


(cos  8  • )  = 


r  (q)-ro(q) 


*  1*2  (q)  ~z5  WUZ+  ir2  {q)  **r5  1 


2  r<<B 
0,  < 


Figure  A. 4,  Definition  of  Normal  to  Solid  Surface 


This  completes  the  definition  of  all  geometrical 
quantities  required  for  the  numerical  treatment  of  par¬ 
tial  cells.  The  scheme  is  used  in  computer  program  SAMS 
to  define  the  partial  cells  associated  with  the  fixed  muz 
zle  boundaries  and  the  moving  projectile  boundaries. 


,  ,  ht  - -  ---  -  " — ■■TS  -.TFP-gng- - —  -,5-  - 


■  q+1 


1st  in-  2nd  index 
dex  of  of  bead 
bead  q  q+1 


j (q+1) = j (q) 


j  (q+1)  >j  (q) 


i (q+1) >i (q) 


If  z4^o  i=; 


If  z4<0 


If  z4sAz(q) 


j (q+1)  j(q) 


If  z4>  Az (q) 


j  (q+) )  “j  (q) 


j (q+1) >j (q) 


i (q+1) <i(q) 


If  z^O 


If  z2<0 


j  (q+1)  <j  (q) 


If  z2s:Az(q) 


If  z2>Az(q) 


i (q+1) =i (q) 


j  (q+1)  >j  (q) 


j  (q+D  <  j  (q) 


Table  A.l.  Definition  of  Boundary  Intersections  With  Cell  Sides 


Type  1-2 

j (q+1) >j (q) 

j (q+1) sj (q) 

Type  1-3 

i 

j  (q+D  >  j  (q) 

j (q+1) sj (q) 

Type  1-4 

i 

f 

j (q+1) > j (q) 

j (q+1) £ j (q) 

Type  2-3 

i (q+1) <i (q) 

i (q+1) si (q) 

Type  2-4 

i (q+l)< i (q) 

1 

| 

i  (q+l)ai  (q) 

Type  3-4 

i (q+1) >i(q) 

i(q+l)si(q) 

r3(q)=r(i(q)-l) 

z3(q)=z(j(q)) 

r4(q)=r3<q) 

z4(q)=z3(q) 


r3  (q)  =r  (i  (q)  -1)  ,  z3  (q)  =2  ( j  (q)  ) 
r4 tq*  *°r3  (q>  <  z4  (<J) =z  ( j  (q)  -i) 


*3  (q)  »r  U  (q)  ) ,  Zj  (q)  =2  ( j  (q)  -1) 
*4  (q)  =r3  (q>  ,  z4  (q)  =z  ( j  (q) ) 


r3(q)=r(i(q) ) 
z3(q)=z(j(q)) 

r4 (q) =r3 (q) 

z4 (q) =z3  Cq) 


r3 (q) =r (i (q) -1) 
z3(q)=z(j (q) -1) 

r4(q)=r3(q) 

z4(q)=z4(q) 


r3(q)“r(i(q)-l)  ,z3(q)=z(j(q)-l 
r4(q)=r3(i(q))  ,z4(q)az,(q) 


r3  (q)  =r  (i  (q)  -X) ,  z3  (q)  »z  ( j  (q) ) 
r4  (q)  Br(i(q)) ,  z4  (q)  «z3  (q) 


r3(q)=r(i(q>) 
z  3(q)=z( j (q)-l) 
r4(q)=r3(q) 
z4(q)=z3(q) 


Table  A. 2.  Final  Definition  of  Partial  Cell  Geometry 


appendix  b 


M-16  MUZZLE  PLOW  TIME  HISTORIES 


The  present  version  of  the  SAMS  code  requires  the 
specification  of  the  time  dependent  muzzle  flow  proper¬ 
ties  (p,  p,  v)  as  an  inflow  boundary  condition  at  some 
location  in  the  barrel  "upstream"  of  the  muzzle  device, 
in  the  calculations  presented  above  for  the  M-16  rifle 
these  were  determined  by  curvefits  to  the  calculated  re¬ 
sults  of  Dr.  Celmins  of  BRL1  .  The  curvefits  are  accu¬ 
rate  to  ±5  percent  for  0  =s  t  £  100  microseconds  and  are 
given  by: 


V (m/sec } 

=  1210  -  110  x  JL 

(B.l) 

P (N/m2 )  = 

(486  -  48.8  x  JL)  x  io5 

(B.2) 

* 

p (kg/m3) 

=  40.95  +  7  x 

/  J  •  J 

(B.  3) 

where  t  is  given  in  microseconds. 


